- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DOM/domvvl.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DOM/domvvl.F90
r12178 r12928 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) … … 208 235 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 236 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / r dt237 frq_rst_hdv(:,:) = 1._wp / rn_Dt 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 / rn_Dt 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 / rn_Dt) & 256 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*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 … … 239 264 ij0 = 128 ; ij1 = 135 ; 240 265 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp 241 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / r dt266 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt 242 267 ENDIF 243 268 ENDIF … … 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 294 320 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 295 REAL(wp) :: z 2dt, z_tmin, z_tmax! local scalars321 REAL(wp) :: z_tmin, z_tmax ! local scalars 296 322 LOGICAL :: ll_do_bclinic ! local logical 297 323 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv … … 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(:,:) ) … … 347 373 IF( kt > nit000 ) THEN 348 374 DO jk = 1, jpkm1 349 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - r dt * frq_rst_hdv(:,:) &350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )375 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 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) … … 414 430 ! 4 - Time stepping of baroclinic scale factors 415 431 ! --------------------------------------------- 416 ! Leapfrog time stepping417 ! ~~~~~~~~~~~~~~~~~~~~~~418 IF( neuler == 0 .AND. kt == nit000 ) THEN419 z2dt = rdt420 ELSE421 z2dt = 2.0_wp * rdt422 ENDIF423 432 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 424 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)433 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 425 434 426 435 ! Maximum deformation control … … 476 485 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 486 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )487 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 488 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)489 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 490 END DO 482 491 … … 486 495 ! ! ---baroclinic part--------- ! 487 496 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)497 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 498 END DO 490 499 ENDIF … … 501 510 zht(:,:) = 0.0_wp 502 511 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(:,:) ) )512 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 513 END DO 514 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 515 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_tmax516 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 517 ! 509 518 zht(:,:) = 0.0_wp 510 519 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(:,:) ) )520 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 521 END DO 522 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 523 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_tmax524 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 525 ! 517 526 zht(:,:) = 0.0_wp 518 527 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(:,:) ) )528 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 529 END DO 530 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 531 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(:,:) ) )532 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 533 ! 534 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 535 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(:,:) ) )536 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 537 ! 538 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 539 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(:,:) ) )540 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 541 ! 542 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 543 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax544 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 545 END IF 537 546 … … 540 549 ! *********************************** ! 541 550 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )551 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 552 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 553 545 554 ! *********************************** ! … … 547 556 ! *********************************** ! 548 557 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)558 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 559 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 560 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)561 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 562 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 563 END DO 555 564 ! ! Inverse of the local depth 556 565 !!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(:,:) )566 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 567 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 568 ! 560 569 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 572 564 573 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***574 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 575 !!---------------------------------------------------------------------- 576 !! *** ROUTINE dom_vvl_sf_update *** 568 577 !! 569 !! ** Purpose : compute time filter and swap of scale factors578 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 579 !! compute all depths and related variables for next time step 571 580 !! write outputs and restart file 572 581 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 582 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 583 !! - reconstruct scale factor at other grid points (interpolate) 575 584 !! - recompute depths and water height fields 576 585 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step586 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 587 !! - Recompute: 579 588 !! e3(u/v)_b 580 !! e3w _n589 !! e3w(:,:,:,Kmm) 581 590 !! e3(u/v)w_b 582 591 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n592 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 593 !! h(u/v) and h(u/v)r 585 594 !! … … 587 596 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 597 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 598 INTEGER, INTENT( in ) :: kt ! time step 599 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 600 ! 591 601 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 605 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 606 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')607 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 608 ! 599 609 IF( kt == nit000 ) THEN 600 610 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'611 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 612 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 613 ENDIF 604 614 ! … … 607 617 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 608 618 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 609 IF( neuler == 0 .AND. kt == nit000) THEN619 IF( l_1st_euler ) THEN 610 620 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 611 621 ELSE 612 622 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 613 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )623 & + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 614 624 ENDIF 615 625 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 626 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 627 624 628 ! Compute all missing vertical scale factor and depths … … 626 630 ! Horizontal scale factor interpolations 627 631 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also632 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 633 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 634 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )635 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 636 633 637 ! 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' )638 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 639 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 640 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 641 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 642 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 643 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 644 641 645 ! 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 646 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 647 gdepw(:,:,1,Kmm) = 0.0_wp 648 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 649 DO_3D_11_11( 2, jpk ) 650 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 651 ! 1 for jk = mikt 652 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 653 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 654 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 655 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 656 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 657 END_3D 658 658 659 659 ! Local depth and Inverse of the local depth of the water 660 660 ! ------------------------------------------------------- 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) 661 ! 662 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 663 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)664 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 665 END DO 668 666 669 667 ! write restart file 670 668 ! ================== 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_ swp669 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 670 ! 671 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 672 ! 673 END SUBROUTINE dom_vvl_sf_update 676 674 677 675 … … 704 702 ! 705 703 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 704 DO_3D_10_10( 1, jpk ) 705 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 706 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 707 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 708 END_3D 715 709 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 716 710 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 717 711 ! 718 712 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 713 DO_3D_10_10( 1, jpk ) 714 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 715 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 716 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 717 END_3D 728 718 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 729 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 730 720 ! 731 721 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 722 DO_3D_10_10( 1, jpk ) 723 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 724 & * r1_e1e2f(ji,jj) & 725 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 726 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 727 END_3D 742 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 743 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 783 769 784 770 785 SUBROUTINE dom_vvl_rst( kt, cdrw )771 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 772 !!--------------------------------------------------------------------- 787 773 !! *** ROUTINE dom_vvl_rst *** … … 795 781 !! they are set to 0. 796 782 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 783 INTEGER , INTENT(in) :: kt ! ocean time-step 784 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 785 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 786 ! 800 787 INTEGER :: ji, jj, jk … … 806 793 IF( ln_rstart ) THEN !* Read the restart file 807 794 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh n, ldxios = lrxios )795 CALL iom_get( numror, jpdom_autoglo, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 796 ! 810 797 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 813 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 814 801 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 802 ! 815 803 ! ! --------- ! 816 804 ! ! all cases ! 817 805 ! ! --------- ! 806 ! 818 807 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 )808 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 809 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 810 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'811 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 812 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)813 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 814 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 815 END WHERE 827 IF( neuler == 0) THEN828 e3t _b(:,:,:) = e3t_n(:,:,:)816 IF( l_1st_euler ) THEN 817 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 818 ENDIF 830 819 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'820 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 821 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 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(:,:,:)836 neuler = 0822 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 824 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 825 l_1st_euler = .true. 837 826 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 828 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 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(:,:,:)843 neuler = 0829 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 830 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 831 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 832 l_1st_euler = .true. 844 833 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'834 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 835 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 IF(lwp) write(numout,*) ' neuler is forced to 0'836 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 848 837 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &838 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 839 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 840 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 841 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)854 neuler = 0842 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 843 l_1st_euler = .true. 855 844 ENDIF 856 845 ! ! ----------- ! … … 888 877 IF( cn_cfg == 'wad' ) THEN 889 878 ! 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 (:,:,:)879 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 880 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 881 ssh (:,:,Kmm) = ssh(:,:,Kbb) 882 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 883 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 884 ELSE 896 885 ! 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 886 ssh(:,:,Kmm) = -ssh_ref 887 ssh(:,:,Kbb) = -ssh_ref 888 889 DO_2D_11_11 890 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 891 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 892 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 893 ENDIF 894 END_2D 910 895 ENDIF !If test case else 911 896 912 897 ! Adjust vertical metrics for all wad 913 898 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &899 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 900 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 901 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 902 END DO 918 e3t_b(:,:,:) = e3t_n(:,:,:) 919 920 DO ji = 1, jpi 921 DO jj = 1, jpj 922 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 923 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 924 ENDIF 925 END DO 926 END DO 903 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 904 905 DO_2D_11_11 906 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 907 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 908 ENDIF 909 END_2D 927 910 ! 928 911 ELSE … … 930 913 ! Just to read set ssh in fact, called latter once vertical grid 931 914 ! is set up: 932 ! CALL usr_def_istate( gdept_0, tmask, ts b, ub, vb, sshb)915 ! CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 933 916 ! ! 934 917 ! DO jk=1,jpk 935 ! e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &918 ! e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 936 919 ! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 937 920 ! END DO 938 ! e3t _n(:,:,:) = e3t_b(:,:,:)939 ssh n(:,:)=0._wp940 e3t _n(:,:,:)=e3t_0(:,:,:)941 e3t _b(:,:,:)=e3t_0(:,:,:)921 ! e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 922 ssh(:,:,Kmm)=0._wp 923 e3t(:,:,:,Kmm)=e3t_0(:,:,:) 924 e3t(:,:,:,Kbb)=e3t_0(:,:,:) 942 925 ! 943 926 END IF ! end of ll_wd edits … … 957 940 ! ! all cases ! 958 941 ! ! --------- ! 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 )942 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 943 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 961 944 ! ! ----------------------- ! 962 945 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 991 974 !!---------------------------------------------------------------------- 992 975 ! 993 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :994 976 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 995 977 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 978 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 998 979 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) … … 1018 999 WRITE(numout,*) ' rn_rst_e3t = 0.e0' 1019 1000 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1020 WRITE(numout,*) ' rn_lf_cutoff = 1.0/r dt'1001 WRITE(numout,*) ' rn_lf_cutoff = 1.0/rn_Dt' 1021 1002 ELSE 1022 1003 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t … … 1033 1014 ! 1034 1015 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 1016 ! 1037 1017 IF(lwp) THEN ! Print the choice
Note: See TracChangeset
for help on using the changeset viewer.
