Newer
Older
MODULE usrdef_zgr
!!======================================================================
!! *** MODULE usrdef_zgr ***
!!
!! === OVERFLOW case ===
!!
!! user defined : vertical coordinate system of a user configuration
!!======================================================================
!! History : 4.0 ! 2016-08 (G. Madec, S. Flavoni) Original code
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! usr_def_zgr : user defined vertical coordinate system (required)
!! zgr_z1d : reference 1D z-coordinate
!!---------------------------------------------------------------------
USE oce ! ocean variables
USE dom_oce ! ocean space and time domain
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
USE usrdef_nam ! User defined : namelist variables
!
USE in_out_manager ! I/O manager
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! distributed memory computing library
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC usr_def_zgr ! called by domzgr.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: usrdef_zgr.F90 15033 2021-06-21 10:24:45Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE usr_def_zgr( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate
& pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate
& pdept , pdepw , & ! 3D t & w-points depth
& pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors
& pe3w , pe3uw , pe3vw, & ! - - -
& k_top , k_bot ) ! top & bottom ocean level
!!---------------------------------------------------------------------
!! *** ROUTINE usr_def_zgr ***
!!
!! ** Purpose : User defined the vertical coordinates
!!
!!----------------------------------------------------------------------
LOGICAL , INTENT(in ) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags ( read in namusr_def )
LOGICAL , INTENT( out) :: ld_isfcav ! under iceshelf cavity flag
REAL(wp), DIMENSION(:) , INTENT( out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:) , INTENT( out) :: pe3t_1d , pe3w_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pdept, pdepw ! grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pe3w , pe3uw, pe3vw ! i-scale factors
INTEGER , DIMENSION(:,:) , INTENT( out) :: k_top, k_bot ! first & last ocean level
!
INTEGER :: ji, jj, jk ! dummy indices
INTEGER :: ik ! local integers
REAL(wp) :: zfact, z1_jpkm1 ! local scalar
REAL(wp) :: ze3min ! local scalar
REAL(wp), DIMENSION(jpi,jpj) :: zht, zhu, z2d ! 2D workspace
!!----------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'usr_def_zgr : OVERFLOW configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
!
!
! type of vertical coordinate
! ---------------------------
! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
ld_isfcav = .FALSE.
!
!
! Build the vertical coordinate system
! ------------------------------------
!
! !== UNmasked meter bathymetry ==!
!
! western continental shelf (500m deep) and eastern deep ocean (2000m deep)
! (set through the jpk and jpi (see userdef_nam.F90))
! with a hyperbolic tangent transition centered at 40km
! NB: here glamt is in kilometers
!
zht(:,:) = + ( 500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,:) - 40.) / 7. ) ) )
!
! at u-point: averaging zht
DO_2D( 0, 0, 0, 0 )
zhu(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji+1,jj) )
END_2D
CALL lbc_lnk( 'usrdef_zgr', zhu, 'U', 1. ) ! boundary condition: this mask the surrouding grid-points
! ! ==>>> set by hand non-zero value on first/last columns & rows
DO ji = mi0(1,nn_hls), mi1(1,nn_hls) ! first row of global domain only
DO ji = mi0(jpiglo,nn_hls), mi1(jpiglo,nn_hls) ! last row of global domain only
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
zhu(ji,2) = zht(ji,2)
END DO
zhu(:,1) = zhu(:,2)
zhu(:,3) = zhu(:,2)
!
CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d ) ! Reference z-coordinate system
!
!
! !== top masked level bathymetry ==! (all coordinates)
!
! no ocean cavities : top ocean level is ONE, except over land
! the ocean basin surrounded by land (1+nn_hls grid-points) set through lbc_lnk call
z2d(:,:) = 1._wp ! surface ocean is the 1st level
CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! closed basin, see userdef_nam.F90
k_top(:,:) = NINT( z2d(:,:) )
!
!
!
IF ( ld_sco ) THEN !== s-coordinate ==! (terrain-following coordinate)
!
k_bot(:,:) = jpkm1 * k_top(:,:) !* bottom ocean = jpk-1 (here use k_top as a land mask)
!
! !* terrain-following coordinate with e3.(k)=cst)
! ! OVERFLOW case : identical with j-index (T=V, U=F)
z1_jpkm1 = 1._wp / REAL( jpkm1 , wp)
DO jk = 1, jpk
pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk , wp ) - 0.5_wp )
pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp ) )
pe3t (:,:,jk) = zht(:,:) * z1_jpkm1
pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1
pe3v (:,:,jk) = zht(:,:) * z1_jpkm1
pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1
pe3w (:,:,jk) = zht(:,:) * z1_jpkm1
pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1
pe3vw(:,:,jk) = zht(:,:) * z1_jpkm1
END DO
ENDIF
!
!
IF ( ld_zco ) THEN !== z-coordinate ==! (step-like topography)
!
! !* bottom ocean compute from the depth of grid-points
k_bot(:,:) = jpkm1 * k_top(:,:) ! here use k_top as a land mask
DO jk = 1, jpkm1
WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) ) k_bot(:,:) = jk * k_top(:,:)
END DO
! !* horizontally uniform coordinate (reference z-co everywhere)
DO jk = 1, jpk
pdept(:,:,jk) = pdept_1d(jk)
pdepw(:,:,jk) = pdepw_1d(jk)
pe3t (:,:,jk) = pe3t_1d (jk)
pe3u (:,:,jk) = pe3t_1d (jk)
pe3v (:,:,jk) = pe3t_1d (jk)
pe3f (:,:,jk) = pe3t_1d (jk)
pe3w (:,:,jk) = pe3w_1d (jk)
pe3uw(:,:,jk) = pe3w_1d (jk)
pe3vw(:,:,jk) = pe3w_1d (jk)
END DO
ENDIF
!
!
IF ( ld_zps ) THEN !== zps-coordinate ==! (partial bottom-steps)
!
ze3min = 0.1_wp * rn_dz
IF(lwp) WRITE(numout,*) ' minimum thickness of the partial cells = 10 % of e3 = ', ze3min
!
!
! !* bottom ocean compute from the depth of grid-points
k_bot(:,:) = jpkm1
DO jk = jpkm1, 1, -1
WHERE( zht(:,:) < pdepw_1d(jk) + ze3min ) k_bot(:,:) = jk-1
END DO
!
! !* vertical coordinate system
DO jk = 1, jpk ! initialization to the reference z-coordinate
pdept(:,:,jk) = pdept_1d(jk)
pdepw(:,:,jk) = pdepw_1d(jk)
pe3t (:,:,jk) = pe3t_1d (jk)
pe3u (:,:,jk) = pe3t_1d (jk)
pe3v (:,:,jk) = pe3t_1d (jk)
pe3f (:,:,jk) = pe3t_1d (jk)
pe3w (:,:,jk) = pe3w_1d (jk)
pe3uw(:,:,jk) = pe3w_1d (jk)
pe3vw(:,:,jk) = pe3w_1d (jk)
END DO
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ik = k_bot(ji,jj)
pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
pe3t (ji,jj,ik ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik )
!
pdept(ji,jj,ik ) = pdepw(ji,jj,ik ) + pe3t (ji,jj,ik ) * 0.5_wp
pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik) ! = pe3t (ji,jj,ik )
pe3w (ji,jj,ik ) = pdept(ji,jj,ik ) - pdept(ji,jj,ik-1) ! st caution ik > 1
END_2D
! ! bottom scale factors and depth at U-, V-, UW and VW-points
! ! usually Computed as the minimum of neighbooring scale factors
pe3u (:,:,:) = pe3t(:,:,:) ! HERE OVERFLOW configuration :
pe3v (:,:,:) = pe3t(:,:,:) ! e3 increases with i-index and identical with j-index
pe3f (:,:,:) = pe3t(:,:,:) ! so e3 minimum of (i,i+1) points is (i) point
pe3uw(:,:,:) = pe3w(:,:,:) ! in j-direction e3v=e3t and e3f=e3v
pe3vw(:,:,:) = pe3w(:,:,:) ! ==>> no need of lbc_lnk calls
!
ENDIF
!
END SUBROUTINE usr_def_zgr
SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d ) ! 1D reference vertical coordinate
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_z1d ***
!!
!! ** Purpose : set the depth of model levels and the resulting
!! vertical scale factors.
!!
!! ** Method : 1D z-coordinate system (use in all type of coordinate)
!! The depth of model levels is set from dep(k), an analytical function:
!! w-level: depw_1d = dep(k)
!! t-level: dept_1d = dep(k+0.5)
!! The scale factors are the discrete derivative of the depth:
!! e3w_1d(jk) = dk[ dept_1d ]
!! e3t_1d(jk) = dk[ depw_1d ]
!!
!! === Here constant vertical resolution ===
!!
!! ** Action : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
!! - pe3t_1d , pe3w_1d : scale factors at T- and W-levels (m)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:), INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:), INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]
!
INTEGER :: jk ! dummy loop indices
REAL(wp) :: zt, zw ! local scalar
!!----------------------------------------------------------------------
!
IF(lwp) THEN ! Parameter print
WRITE(numout,*)
WRITE(numout,*) ' zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
WRITE(numout,*) ' ~~~~~~~'
ENDIF
!
! Reference z-coordinate (depth - scale factor at T- and W-points) ! Madec & Imbard 1996 function
! ----------------------
DO jk = 1, jpk
zw = REAL( jk , wp )
zt = REAL( jk , wp ) + 0.5_wp
pdepw_1d(jk) = rn_dz * REAL( jk-1 , wp )
pdept_1d(jk) = rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
pe3w_1d (jk) = rn_dz
pe3t_1d (jk) = rn_dz
END DO
!
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:'
WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
ENDIF
!
END SUBROUTINE zgr_z1d
!!======================================================================
END MODULE usrdef_zgr