Changeset 12377 for NEMO/trunk/src/OCE/DOM/domvvl.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DOM/domvvl.F90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/DOM/domvvl.F90
r11536 r12377 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 10 11 !!---------------------------------------------------------------------- 11 12 … … 13 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 14 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 !! dom_vvl_sf_ swp: Swap vertical scale factors and update the vertical grid16 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 16 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 17 18 !! dom_vvl_rst : read/write restart file … … 36 37 37 38 PUBLIC dom_vvl_init ! called by domain.F90 39 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 38 40 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_ swp! called by step.F9041 PUBLIC dom_vvl_sf_update ! called by step.F90 40 42 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 43 … … 62 64 63 65 !! * Substitutions 64 # include " vectopt_loop_substitute.h90"66 # include "do_loop_substitute.h90" 65 67 !!---------------------------------------------------------------------- 66 68 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 93 95 94 96 95 SUBROUTINE dom_vvl_init 97 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 98 !!---------------------------------------------------------------------- 97 99 !! *** ROUTINE dom_vvl_init *** … … 102 104 !! ** Method : - use restart file and/or initialize 103 105 !! - interpolate scale factors 106 !! 107 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 108 !! - Regrid: e3[u/v](:,:,:,Kmm) 109 !! e3[u/v](:,:,:,Kmm) 110 !! e3w(:,:,:,Kmm) 111 !! e3[u/v]w_b 112 !! e3[u/v]w_n 113 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 114 !! - h(t/u/v)_0 115 !! - frq_rst_e3t and frq_rst_hdv 116 !! 117 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 118 !!---------------------------------------------------------------------- 119 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 120 ! 121 IF(lwp) WRITE(numout,*) 122 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 123 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 124 ! 125 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 126 ! 127 ! ! Allocate module arrays 128 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 129 ! 130 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 131 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 132 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 133 ! 134 CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 135 ! 136 END SUBROUTINE dom_vvl_init 137 ! 138 SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 139 !!---------------------------------------------------------------------- 140 !! *** ROUTINE dom_vvl_init *** 141 !! 142 !! ** Purpose : Interpolation of all scale factors, 143 !! depths and water column heights 144 !! 145 !! ** Method : - interpolate scale factors 104 146 !! 105 147 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) … … 115 157 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 158 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 160 !!---------------------------------------------------------------------- 117 161 INTEGER :: ji, jj, jk 118 162 INTEGER :: ii0, ii1, ij0, ij1 … … 120 164 !!---------------------------------------------------------------------- 121 165 ! 122 IF(lwp) WRITE(numout,*)123 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'124 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'125 !126 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)127 !128 ! ! Allocate module arrays129 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )130 !131 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf132 CALL dom_vvl_rst( nit000, 'READ' )133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all134 !135 166 ! !== Set of all other vertical scale factors ==! (now and before) 136 167 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t _b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U138 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )139 CALL dom_vvl_interpol( e3t _b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V140 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )141 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F168 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 169 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 170 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 171 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 172 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 173 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W144 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b (:,:,:), 'W' )145 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW146 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )147 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW148 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )174 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 175 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 176 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 177 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 178 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 179 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 180 150 181 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t _a(:,:,:) = e3t_n(:,:,:)152 e3u _a(:,:,:) = e3u_n(:,:,:)153 e3v _a(:,:,:) = e3v_n(:,:,:)182 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 183 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 184 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 185 ! 155 186 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 156 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 157 gdepw_n(:,:,1) = 0.0_wp 158 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 159 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 160 gdepw_b(:,:,1) = 0.0_wp 161 DO jk = 2, jpk ! vertical sum 162 DO jj = 1,jpj 163 DO ji = 1,jpi 164 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 165 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 168 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 169 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 170 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 171 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 172 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 173 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 174 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 175 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 176 END DO 177 END DO 187 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 188 gdepw(:,:,1,Kmm) = 0.0_wp 189 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 190 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 191 gdepw(:,:,1,Kbb) = 0.0_wp 192 DO_3D_11_11( 2, jpk ) 193 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 194 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 195 ! ! 0.5 where jk = mikt 196 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 197 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 198 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 199 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 200 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 201 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 202 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 203 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 204 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 205 END_3D 206 ! 207 ! !== thickness of the water column !! (ocean portion only) 208 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 209 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 210 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 211 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 212 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 213 DO jk = 2, jpkm1 214 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 215 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 216 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 217 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 218 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 178 219 END DO 179 220 ! 180 ! !== thickness of the water column !! (ocean portion only)181 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....182 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)183 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)184 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)185 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)186 DO jk = 2, jpkm1187 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)188 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)189 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)190 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)191 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)192 END DO193 !194 221 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF196 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )197 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )198 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )222 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 223 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 224 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 225 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 199 226 200 227 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 211 238 ENDIF 212 239 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 213 DO jj = 1, jpj 214 DO ji = 1, jpi 240 DO_2D_11_11 215 241 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 216 IF( ABS(gphit(ji,jj)) >= 6.) THEN 217 ! values outside the equatorial band and transition zone (ztilde) 218 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 219 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 220 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 221 ! values inside the equatorial band (ztilde as zstar) 222 frq_rst_e3t(ji,jj) = 0.0_wp 223 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 224 ELSE ! transition band (2.5 to 6 degrees N/S) 225 ! ! (linearly transition from z-tilde to z-star) 226 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 227 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 228 & * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 230 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 231 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 232 & * 180._wp / 3.5_wp ) ) 233 ENDIF 234 END DO 235 END DO 242 IF( ABS(gphit(ji,jj)) >= 6.) THEN 243 ! values outside the equatorial band and transition zone (ztilde) 244 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 245 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 246 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 247 ! values inside the equatorial band (ztilde as zstar) 248 frq_rst_e3t(ji,jj) = 0.0_wp 249 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 250 ELSE ! transition band (2.5 to 6 degrees N/S) 251 ! ! (linearly transition from z-tilde to z-star) 252 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 253 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 254 & * 180._wp / 3.5_wp ) ) 255 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 256 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 257 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 258 & * 180._wp / 3.5_wp ) ) 259 ENDIF 260 END_2D 236 261 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 237 262 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 … … 263 288 ENDIF 264 289 ! 265 END SUBROUTINE dom_vvl_ init266 267 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )290 END SUBROUTINE dom_vvl_zgr 291 292 293 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 294 !!---------------------------------------------------------------------- 270 295 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 313 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 314 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 315 INTEGER, INTENT( in ) :: kt ! time step 316 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 317 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 318 ! 293 319 INTEGER :: ji, jj, jk ! dummy loop indices … … 321 347 ! ! --------------------------------------------- ! 322 348 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )349 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 350 DO jk = 1, jpkm1 325 ! formally this is the same as e3t _a= e3t_0*(1+ssha/ht_0)326 e3t _a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)351 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 352 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 353 END DO 328 354 ! … … 337 363 zht(:,:) = 0._wp 338 364 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)365 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 366 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 367 END DO 342 368 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 348 374 DO jk = 1, jpkm1 349 375 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )376 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 377 END DO 352 378 ENDIF … … 361 387 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 388 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )389 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 390 END DO 365 391 ELSE ! layer case 366 392 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)393 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 394 END DO 369 395 ENDIF … … 381 407 zwu(:,:) = 0._wp 382 408 zwv(:,:) = 0._wp 383 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 384 DO jj = 1, jpjm1 385 DO ji = 1, fs_jpim1 ! vector opt. 386 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 387 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 388 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 389 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 390 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 391 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 392 END DO 393 END DO 394 END DO 395 DO jj = 1, jpj ! b - correction for last oceanic u-v points 396 DO ji = 1, jpi 397 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 398 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 399 END DO 400 END DO 401 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vector opt. 404 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 405 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 406 & ) * r1_e1e2t(ji,jj) 407 END DO 408 END DO 409 END DO 409 DO_3D_10_10( 1, jpkm1 ) 410 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 411 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 412 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 413 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 414 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 415 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 416 END_3D 417 DO_2D_11_11 418 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 419 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 420 END_2D 421 DO_3D_00_00( 1, jpkm1 ) 422 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 423 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 424 & ) * r1_e1e2t(ji,jj) 425 END_3D 410 426 ! ! d - thickness diffusion transport: boundary conditions 411 427 ! (stored for tracer advction and continuity equation) … … 476 492 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 493 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )494 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 495 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)496 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 497 END DO 482 498 … … 486 502 ! ! ---baroclinic part--------- ! 487 503 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)504 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 505 END DO 490 506 ENDIF … … 501 517 zht(:,:) = 0.0_wp 502 518 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)504 END DO 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh n(:,:) - zht(:,:) ) )519 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 520 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 522 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t _n))) =', z_tmax523 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 524 ! 509 525 zht(:,:) = 0.0_wp 510 526 DO jk = 1, jpkm1 511 zht(:,:) = zht(:,:) + e3t _a(:,:,jk) * tmask(:,:,jk)512 END DO 513 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh a(:,:) - zht(:,:) ) )527 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 528 END DO 529 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 530 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 515 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t _a))) =', z_tmax531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 532 ! 517 533 zht(:,:) = 0.0_wp 518 534 DO jk = 1, jpkm1 519 zht(:,:) = zht(:,:) + e3t _b(:,:,jk) * tmask(:,:,jk)520 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh b(:,:) - zht(:,:) ) )535 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 536 END DO 537 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 538 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 523 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t _b))) =', z_tmax524 ! 525 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh b(:,:) ) )539 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 540 ! 541 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 542 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh b))) =', z_tmax528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh n(:,:) ) )543 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 544 ! 545 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 546 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh n))) =', z_tmax532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh a(:,:) ) )547 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 548 ! 549 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 550 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax551 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 552 END IF 537 553 … … 540 556 ! *********************************** ! 541 557 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )558 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 559 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 560 545 561 ! *********************************** ! … … 547 563 ! *********************************** ! 548 564 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)565 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 566 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 567 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)568 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 569 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 570 END DO 555 571 ! ! Inverse of the local depth 556 572 !!gm BUG ? don't understand the use of umask_i here ..... 557 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )558 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )573 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 574 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 575 ! 560 576 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 579 564 580 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***581 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 582 !!---------------------------------------------------------------------- 583 !! *** ROUTINE dom_vvl_sf_update *** 568 584 !! 569 !! ** Purpose : compute time filter and swap of scale factors585 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 586 !! compute all depths and related variables for next time step 571 587 !! write outputs and restart file 572 588 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 589 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 590 !! - reconstruct scale factor at other grid points (interpolate) 575 591 !! - recompute depths and water height fields 576 592 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step593 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 594 !! - Recompute: 579 595 !! e3(u/v)_b 580 !! e3w _n596 !! e3w(:,:,:,Kmm) 581 597 !! e3(u/v)w_b 582 598 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n599 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 600 !! h(u/v) and h(u/v)r 585 601 !! … … 587 603 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 604 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 605 INTEGER, INTENT( in ) :: kt ! time step 606 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 607 ! 591 608 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 612 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 613 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')614 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 615 ! 599 616 IF( kt == nit000 ) THEN 600 617 IF(lwp) WRITE(numout,*) 601 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'602 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'618 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 619 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 620 ENDIF 604 621 ! … … 615 632 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 633 ENDIF 617 gdept_b(:,:,:) = gdept_n(:,:,:)618 gdepw_b(:,:,:) = gdepw_n(:,:,:)619 620 e3t_n(:,:,:) = e3t_a(:,:,:)621 e3u_n(:,:,:) = e3u_a(:,:,:)622 e3v_n(:,:,:) = e3v_a(:,:,:)623 634 624 635 ! Compute all missing vertical scale factor and depths … … 626 637 ! Horizontal scale factor interpolations 627 638 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also639 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 640 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 641 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )642 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 643 633 644 ! Vertical scale factor interpolations 634 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n(:,:,:), 'W' )635 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' )636 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' )637 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b(:,:,:), 'W' )638 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )639 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )645 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 646 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 647 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 648 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 649 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 650 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 651 641 652 ! t- and w- points depth (set the isf depth as it is in the initial step) 642 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 643 gdepw_n(:,:,1) = 0.0_wp 644 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 645 DO jk = 2, jpk 646 DO jj = 1,jpj 647 DO ji = 1,jpi 648 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 649 ! 1 for jk = mikt 650 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 651 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 652 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 653 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 654 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 655 END DO 656 END DO 657 END DO 653 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 654 gdepw(:,:,1,Kmm) = 0.0_wp 655 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 656 DO_3D_11_11( 2, jpk ) 657 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 658 ! 1 for jk = mikt 659 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 660 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 661 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 662 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 663 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 664 END_3D 658 665 659 666 ! Local depth and Inverse of the local depth of the water 660 667 ! ------------------------------------------------------- 661 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 662 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 ! 664 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 668 ! 669 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 670 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)671 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 672 END DO 668 673 669 674 ! write restart file 670 675 ! ================== 671 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )672 ! 673 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')674 ! 675 END SUBROUTINE dom_vvl_sf_ swp676 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 677 ! 678 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 679 ! 680 END SUBROUTINE dom_vvl_sf_update 676 681 677 682 … … 704 709 ! 705 710 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 706 DO jk = 1, jpk 707 DO jj = 1, jpjm1 708 DO ji = 1, fs_jpim1 ! vector opt. 709 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 710 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 711 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 712 END DO 713 END DO 714 END DO 711 DO_3D_10_10( 1, jpk ) 712 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 713 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 714 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 715 END_3D 715 716 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 716 717 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 717 718 ! 718 719 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 719 DO jk = 1, jpk 720 DO jj = 1, jpjm1 721 DO ji = 1, fs_jpim1 ! vector opt. 722 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 723 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 724 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 725 END DO 726 END DO 727 END DO 720 DO_3D_10_10( 1, jpk ) 721 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 722 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 723 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 724 END_3D 728 725 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 729 726 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 730 727 ! 731 728 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 732 DO jk = 1, jpk 733 DO jj = 1, jpjm1 734 DO ji = 1, fs_jpim1 ! vector opt. 735 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 736 & * r1_e1e2f(ji,jj) & 737 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 738 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 739 END DO 740 END DO 741 END DO 729 DO_3D_10_10( 1, jpk ) 730 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 731 & * r1_e1e2f(ji,jj) & 732 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 733 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 734 END_3D 742 735 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 743 736 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 783 776 784 777 785 SUBROUTINE dom_vvl_rst( kt, cdrw )778 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 779 !!--------------------------------------------------------------------- 787 780 !! *** ROUTINE dom_vvl_rst *** … … 795 788 !! they are set to 0. 796 789 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 790 INTEGER , INTENT(in) :: kt ! ocean time-step 791 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 792 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 793 ! 800 794 INTEGER :: ji, jj, jk … … 806 800 IF( ln_rstart ) THEN !* Read the restart file 807 801 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh n, ldxios = lrxios )802 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 803 ! 810 804 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 813 807 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 814 808 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 809 ! 815 810 ! ! --------- ! 816 811 ! ! all cases ! 817 812 ! ! --------- ! 813 ! 818 814 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 819 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )820 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )815 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 816 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 817 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'818 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 819 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)820 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 821 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 822 END WHERE 827 823 IF( neuler == 0 ) THEN 828 e3t _b(:,:,:) = e3t_n(:,:,:)824 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 825 ENDIF 830 826 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 828 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 829 IF(lwp) write(numout,*) 'neuler is forced to 0' 834 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t _b(:,:,:), ldxios = lrxios )835 e3t _n(:,:,:) = e3t_b(:,:,:)830 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 831 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 836 832 neuler = 0 837 833 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'834 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 835 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 836 IF(lwp) write(numout,*) 'neuler is forced to 0' 841 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t _n(:,:,:), ldxios = lrxios )842 e3t _b(:,:,:) = e3t_n(:,:,:)837 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 838 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 843 839 neuler = 0 844 840 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'841 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 842 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 843 IF(lwp) write(numout,*) 'neuler is forced to 0' 848 844 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &845 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 846 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 847 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 848 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)849 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 854 850 neuler = 0 855 851 ENDIF … … 888 884 IF( cn_cfg == 'wad' ) THEN 889 885 ! Wetting and drying test case 890 CALL usr_def_istate( gdept _b, tmask, tsb, ub, vb, sshb)891 ts n (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones892 ssh n (:,:) = sshb(:,:)893 u n (:,:,:) = ub (:,:,:)894 v n (:,:,:) = vb (:,:,:)886 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 887 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 888 ssh (:,:,Kmm) = ssh(:,:,Kbb) 889 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 890 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 891 ELSE 896 892 ! if not test case 897 sshn(:,:) = -ssh_ref 898 sshb(:,:) = -ssh_ref 899 900 DO jj = 1, jpj 901 DO ji = 1, jpi 902 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 903 904 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 906 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 907 ENDIF 908 ENDDO 909 ENDDO 893 ssh(:,:,Kmm) = -ssh_ref 894 ssh(:,:,Kbb) = -ssh_ref 895 896 DO_2D_11_11 897 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 898 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 899 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 900 ENDIF 901 END_2D 910 902 ENDIF !If test case else 911 903 912 904 ! Adjust vertical metrics for all wad 913 905 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &906 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 907 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 908 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 909 END DO 918 e3t _b(:,:,:) = e3t_n(:,:,:)910 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 919 911 920 912 DO ji = 1, jpi … … 930 922 ! Just to read set ssh in fact, called latter once vertical grid 931 923 ! is set up: 932 ! CALL usr_def_istate( gdept_0, tmask, ts b, ub, vb, sshb)924 ! CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 933 925 ! ! 934 926 ! DO jk=1,jpk 935 ! e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &927 ! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 936 928 ! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 937 929 ! END DO 938 ! e3t _n(:,:,:) = e3t_b(:,:,:)939 ssh n(:,:)=0._wp940 e3t _n(:,:,:)=e3t_0(:,:,:)941 e3t _b(:,:,:)=e3t_0(:,:,:)930 ! e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 931 ssh(:,:,Kmm)=0._wp 932 e3t(:,:,:,Kmm)=e3t_0(:,:,:) 933 e3t(:,:,:,Kbb)=e3t_0(:,:,:) 942 934 ! 943 935 END IF ! end of ll_wd edits … … 957 949 ! ! all cases ! 958 950 ! ! --------- ! 959 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t _b(:,:,:), ldxios = lwxios )960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t _n(:,:,:), ldxios = lwxios )951 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 952 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 961 953 ! ! ----------------------- ! 962 954 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 991 983 !!---------------------------------------------------------------------- 992 984 ! 993 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :994 985 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 995 986 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' ) 996 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run997 987 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 998 988 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) … … 1033 1023 ! 1034 1024 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1035 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1036 1025 ! 1037 1026 IF(lwp) THEN ! Print the choice
Note: See TracChangeset
for help on using the changeset viewer.
