- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
tests/CANAL/MY_SRC/domvvl.F90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/CANAL/MY_SRC/domvvl.F90
r10425 r13463 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 11 !! 4.x ! 2020-02 (G. Madec, S. Techene) introduce ssh to h0 ratio 10 12 !!---------------------------------------------------------------------- 11 13 12 !!----------------------------------------------------------------------13 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness14 !! dom_vvl_sf_nxt : Compute next vertical scale factors15 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another17 !! dom_vvl_rst : read/write restart file18 !! dom_vvl_ctl : Check the vvl options19 !!----------------------------------------------------------------------20 14 USE oce ! ocean dynamics and tracers 21 15 USE phycst ! physical constant … … 35 29 PRIVATE 36 30 37 PUBLIC dom_vvl_init ! called by domain.F9038 PUBLIC dom_vvl_sf_nxt ! called by step.F9039 PUBLIC dom_vvl_sf_swp ! called by step.F9040 PUBLIC dom_vvl_interpol ! called by dynnxt.F9041 42 31 ! !!* Namelist nam_vvl 43 32 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate … … 61 50 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 62 51 52 #if defined key_qco 53 !!---------------------------------------------------------------------- 54 !! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate 55 !!---------------------------------------------------------------------- 56 #else 57 !!---------------------------------------------------------------------- 58 !! Default key Old management of time varying vertical coordinate 59 !!---------------------------------------------------------------------- 60 61 !!---------------------------------------------------------------------- 62 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 63 !! dom_vvl_sf_nxt : Compute next vertical scale factors 64 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 65 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 66 !! dom_vvl_rst : read/write restart file 67 !! dom_vvl_ctl : Check the vvl options 68 !!---------------------------------------------------------------------- 69 70 PUBLIC dom_vvl_init ! called by domain.F90 71 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 72 PUBLIC dom_vvl_sf_nxt ! called by step.F90 73 PUBLIC dom_vvl_sf_update ! called by step.F90 74 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 75 63 76 !! * Substitutions 64 # include " vectopt_loop_substitute.h90"77 # include "do_loop_substitute.h90" 65 78 !!---------------------------------------------------------------------- 66 79 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 93 106 94 107 95 SUBROUTINE dom_vvl_init 108 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 109 !!---------------------------------------------------------------------- 97 110 !! *** ROUTINE dom_vvl_init *** … … 102 115 !! ** Method : - use restart file and/or initialize 103 116 !! - interpolate scale factors 117 !! 118 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 119 !! - Regrid: e3[u/v](:,:,:,Kmm) 120 !! e3[u/v](:,:,:,Kmm) 121 !! e3w(:,:,:,Kmm) 122 !! e3[u/v]w_b 123 !! e3[u/v]w_n 124 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 125 !! - h(t/u/v)_0 126 !! - frq_rst_e3t and frq_rst_hdv 127 !! 128 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 129 !!---------------------------------------------------------------------- 130 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 131 ! 132 IF(lwp) WRITE(numout,*) 133 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 134 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 135 ! 136 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 137 ! 138 ! ! Allocate module arrays 139 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 140 ! 141 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 142 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 143 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 144 ! 145 CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 146 ! 147 END SUBROUTINE dom_vvl_init 148 149 150 SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 151 !!---------------------------------------------------------------------- 152 !! *** ROUTINE dom_vvl_init *** 153 !! 154 !! ** Purpose : Interpolation of all scale factors, 155 !! depths and water column heights 156 !! 157 !! ** Method : - interpolate scale factors 104 158 !! 105 159 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) … … 115 169 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 170 !!---------------------------------------------------------------------- 171 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 172 !!---------------------------------------------------------------------- 117 173 INTEGER :: ji, jj, jk 118 174 INTEGER :: ii0, ii1, ij0, ij1 … … 120 176 !!---------------------------------------------------------------------- 121 177 ! 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 178 ! !== Set of all other vertical scale factors ==! (now and before) 136 179 ! ! 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 F180 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 181 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 182 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 183 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 184 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 185 ! ! 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' )186 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 187 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 188 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 189 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 190 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 191 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 192 150 193 ! 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(:,:,:)194 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 195 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 196 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 197 ! 155 198 ! !== 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 199 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 200 gdepw(:,:,1,Kmm) = 0.0_wp 201 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 202 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 203 gdepw(:,:,1,Kbb) = 0.0_wp 204 DO_3D( 1, 1, 1, 1, 2, jpk ) 205 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 206 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 207 ! ! 0.5 where jk = mikt 208 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 209 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 210 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 211 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 212 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 213 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 214 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 215 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 216 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 217 END_3D 218 ! 219 ! !== thickness of the water column !! (ocean portion only) 220 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 221 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 222 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 223 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 224 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 225 DO jk = 2, jpkm1 226 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 227 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 228 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 229 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 230 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 178 231 END DO 179 232 ! 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 233 ! !== 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(:,:) )234 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 235 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 236 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 237 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 199 238 200 239 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 208 247 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 248 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / r dt249 frq_rst_hdv(:,:) = 1._wp / rn_Dt 211 250 ENDIF 212 251 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 252 DO_2D( 1, 1, 1, 1 ) 215 253 !!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 254 IF( ABS(gphit(ji,jj)) >= 6.) THEN 255 ! values outside the equatorial band and transition zone (ztilde) 256 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 257 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 258 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 259 ! values inside the equatorial band (ztilde as zstar) 260 frq_rst_e3t(ji,jj) = 0.0_wp 261 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 262 ELSE ! transition band (2.5 to 6 degrees N/S) 263 ! ! (linearly transition from z-tilde to z-star) 264 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 265 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 266 & * 180._wp / 3.5_wp ) ) 267 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 268 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 269 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 270 & * 180._wp / 3.5_wp ) ) 271 ENDIF 272 END_2D 236 273 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 237 274 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 238 ii0 = 103 ; ii1 = 111239 ij0 = 128 ; ij1 = 135 ;275 ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1 276 ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls 240 277 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 dt278 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt 242 279 ENDIF 243 280 ENDIF … … 263 300 ENDIF 264 301 ! 265 END SUBROUTINE dom_vvl_ init266 267 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )302 END SUBROUTINE dom_vvl_zgr 303 304 305 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 306 !!---------------------------------------------------------------------- 270 307 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 325 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 326 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 327 INTEGER, INTENT( in ) :: kt ! time step 328 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 329 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 330 ! 293 331 INTEGER :: ji, jj, jk ! dummy loop indices 294 332 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 295 REAL(wp) :: z 2dt, z_tmin, z_tmax! local scalars333 REAL(wp) :: z_tmin, z_tmax ! local scalars 296 334 LOGICAL :: ll_do_bclinic ! local logical 297 335 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 298 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 336 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3t 337 LOGICAL , DIMENSION(:,:,:), ALLOCATABLE :: llmsk 299 338 !!---------------------------------------------------------------------- 300 339 ! … … 321 360 ! ! --------------------------------------------- ! 322 361 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )362 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 363 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)364 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 365 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 366 END DO 328 367 ! 329 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !330 ! ! ------baroclinic part------ !368 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 369 ! ! ------baroclinic part------ ! 331 370 ! I - initialization 332 371 ! ================== … … 337 376 zht(:,:) = 0._wp 338 377 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)378 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 379 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 380 END DO 342 381 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 347 386 IF( kt > nit000 ) THEN 348 387 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(:,:) ) )388 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 389 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 390 END DO 352 391 ENDIF … … 361 400 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 401 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )402 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 403 END DO 365 404 ELSE ! layer case 366 405 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)406 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 407 END DO 369 408 ENDIF … … 381 420 zwu(:,:) = 0._wp 382 421 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 422 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 423 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 424 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 425 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 426 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 427 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 428 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 429 END_3D 430 DO_2D( 1, 1, 1, 1 ) 431 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 432 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 433 END_2D 434 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 435 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 436 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 437 & ) * r1_e1e2t(ji,jj) 438 END_3D 410 439 ! ! d - thickness diffusion transport: boundary conditions 411 440 ! (stored for tracer advction and continuity equation) … … 414 443 ! 4 - Time stepping of baroclinic scale factors 415 444 ! --------------------------------------------- 416 ! Leapfrog time stepping417 ! ~~~~~~~~~~~~~~~~~~~~~~418 IF( neuler == 0 .AND. kt == nit000 ) THEN419 z2dt = rdt420 ELSE421 z2dt = 2.0_wp * rdt422 ENDIF423 445 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 424 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)446 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 425 447 426 448 ! Maximum deformation control 427 449 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 428 ze3t(:,:,jpk) = 0._wp 429 DO jk = 1, jpkm1 430 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 431 END DO 432 z_tmax = MAXVAL( ze3t(:,:,:) ) 433 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 434 z_tmin = MINVAL( ze3t(:,:,:) ) 435 CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain 450 ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) ) 451 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 452 ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 453 END_3D 454 ! 455 llmsk( 1:Nis1,:,:) = .FALSE. ! exclude halos from the checked region 456 llmsk(Nie1: jpi,:,:) = .FALSE. 457 llmsk(:, 1:Njs1,:) = .FALSE. 458 llmsk(:,Nje1: jpj,:) = .FALSE. 459 ! 460 llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp ! define only the inner domain 461 z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 462 z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain 436 463 ! - ML - test: for the moment, stop simulation for too large e3_t variations 437 464 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 438 IF( lk_mpp ) THEN 439 CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max ) 440 CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min ) 441 ELSE 442 ijk_max = MAXLOC( ze3t(:,:,:) ) 443 ijk_max(1) = ijk_max(1) + nimpp - 1 444 ijk_max(2) = ijk_max(2) + njmpp - 1 445 ijk_min = MINLOC( ze3t(:,:,:) ) 446 ijk_min(1) = ijk_min(1) + nimpp - 1 447 ijk_min(2) = ijk_min(2) + njmpp - 1 448 ENDIF 465 CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max ) 466 CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min ) 449 467 IF (lwp) THEN 450 468 WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax … … 455 473 ENDIF 456 474 ENDIF 475 DEALLOCATE( ze3t, llmsk ) 457 476 ! - ML - end test 458 477 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below … … 476 495 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 496 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )497 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 498 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)499 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 500 END DO 482 501 … … 486 505 ! ! ---baroclinic part--------- ! 487 506 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)507 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 508 END DO 490 509 ENDIF … … 501 520 zht(:,:) = 0.0_wp 502 521 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)522 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 504 523 END DO 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh n(:,:) - zht(:,:) ) )524 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 525 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_tmax526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 527 ! 509 528 zht(:,:) = 0.0_wp 510 529 DO jk = 1, jpkm1 511 zht(:,:) = zht(:,:) + e3t _a(:,:,jk) * tmask(:,:,jk)530 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 512 531 END DO 513 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh a(:,:) - zht(:,:) ) )532 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 533 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_tmax534 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 535 ! 517 536 zht(:,:) = 0.0_wp 518 537 DO jk = 1, jpkm1 519 zht(:,:) = zht(:,:) + e3t _b(:,:,jk) * tmask(:,:,jk)538 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 520 539 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh b(:,:) - zht(:,:) ) )540 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 541 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(:,:) ) )542 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 543 ! 544 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 545 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(:,:) ) )546 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 547 ! 548 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 549 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(:,:) ) )550 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 551 ! 552 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 553 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax554 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 555 END IF 537 556 … … 540 559 ! *********************************** ! 541 560 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )561 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 562 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 563 545 564 ! *********************************** ! … … 547 566 ! *********************************** ! 548 567 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)568 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 569 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 570 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)571 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 572 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 573 END DO 555 574 ! ! Inverse of the local depth 556 575 !!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(:,:) )576 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 577 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 578 ! 560 579 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 582 564 583 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***584 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 585 !!---------------------------------------------------------------------- 586 !! *** ROUTINE dom_vvl_sf_update *** 568 587 !! 569 !! ** Purpose : compute time filter and swap of scale factors588 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 589 !! compute all depths and related variables for next time step 571 590 !! write outputs and restart file 572 591 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 592 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 593 !! - reconstruct scale factor at other grid points (interpolate) 575 594 !! - recompute depths and water height fields 576 595 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step596 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 597 !! - Recompute: 579 598 !! e3(u/v)_b 580 !! e3w _n599 !! e3w(:,:,:,Kmm) 581 600 !! e3(u/v)w_b 582 601 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n602 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 603 !! h(u/v) and h(u/v)r 585 604 !! … … 587 606 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 607 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 608 INTEGER, INTENT( in ) :: kt ! time step 609 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 610 ! 591 611 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 615 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 616 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')617 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 618 ! 599 619 IF( kt == nit000 ) THEN 600 620 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'621 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 622 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 623 ENDIF 604 624 ! … … 607 627 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 608 628 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 609 IF( neuler == 0 .AND. kt == nit000) THEN629 IF( l_1st_euler ) THEN 610 630 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 611 631 ELSE 612 632 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 613 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )633 & + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 614 634 ENDIF 615 635 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 636 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 637 624 638 ! Compute all missing vertical scale factor and depths … … 626 640 ! Horizontal scale factor interpolations 627 641 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also642 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 643 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 644 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )645 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 646 633 647 ! 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' )648 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 649 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 650 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 651 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 652 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 653 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 654 641 655 ! 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 656 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 657 gdepw(:,:,1,Kmm) = 0.0_wp 658 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 659 DO_3D( 1, 1, 1, 1, 2, jpk ) 660 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 661 ! 1 for jk = mikt 662 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 663 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 664 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 665 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 666 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 667 END_3D 658 668 659 669 ! Local depth and Inverse of the local depth of the water 660 670 ! ------------------------------------------------------- 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) 671 ! 672 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 673 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)674 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 675 END DO 668 676 669 677 ! write restart file 670 678 ! ================== 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_ swp679 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 680 ! 681 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 682 ! 683 END SUBROUTINE dom_vvl_sf_update 676 684 677 685 … … 704 712 ! 705 713 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 714 DO_3D( 1, 0, 1, 0, 1, jpk ) 715 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 716 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 717 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 718 END_3D 715 719 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 716 720 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 717 721 ! 718 722 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 723 DO_3D( 1, 0, 1, 0, 1, jpk ) 724 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 725 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 726 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 727 END_3D 728 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 729 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 730 730 ! 731 731 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 732 DO_3D( 1, 0, 1, 0, 1, jpk ) 733 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 734 & * r1_e1e2f(ji,jj) & 735 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 736 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 737 END_3D 742 738 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 743 739 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 783 779 784 780 785 SUBROUTINE dom_vvl_rst( kt, cdrw )781 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 782 !!--------------------------------------------------------------------- 787 783 !! *** ROUTINE dom_vvl_rst *** … … 795 791 !! they are set to 0. 796 792 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 793 INTEGER , INTENT(in) :: kt ! ocean time-step 794 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 795 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 796 ! 800 797 INTEGER :: ji, jj, jk … … 806 803 IF( ln_rstart ) THEN !* Read the restart file 807 804 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_auto glo, 'sshn' , sshn, ldxios = lrxios )805 CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 806 ! 810 807 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 813 810 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 814 811 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 812 ! 815 813 ! ! --------- ! 816 814 ! ! all cases ! 817 815 ! ! --------- ! 816 ! 818 817 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 819 CALL iom_get( numror, jpdom_auto glo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )820 CALL iom_get( numror, jpdom_auto glo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )818 CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 819 CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 820 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'821 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 822 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)823 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 824 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 825 END WHERE 827 IF( neuler == 0) THEN828 e3t _b(:,:,:) = e3t_n(:,:,:)826 IF( l_1st_euler ) THEN 827 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 828 ENDIF 830 829 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'830 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 831 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_auto glo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )835 e3t _n(:,:,:) = e3t_b(:,:,:)836 neuler = 0832 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 833 CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 834 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 835 l_1st_euler = .true. 837 836 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'837 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 838 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_auto glo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )842 e3t _b(:,:,:) = e3t_n(:,:,:)843 neuler = 0839 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 840 CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 841 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 842 l_1st_euler = .true. 844 843 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'844 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 845 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 IF(lwp) write(numout,*) ' neuler is forced to 0'846 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 848 847 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &848 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 849 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 850 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 851 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)854 neuler = 0852 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 853 l_1st_euler = .true. 855 854 ENDIF 856 855 ! ! ----------- ! … … 864 863 ! ! ----------------------- ! 865 864 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 866 CALL iom_get( numror, jpdom_auto glo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )867 CALL iom_get( numror, jpdom_auto glo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios ) 866 CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios ) 868 867 ELSE ! one at least array is missing 869 868 tilde_e3t_b(:,:,:) = 0.0_wp … … 874 873 ! ! ------------ ! 875 874 IF( id5 > 0 ) THEN ! required array exists 876 CALL iom_get( numror, jpdom_auto glo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )875 CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 877 876 ELSE ! array is missing 878 877 hdiv_lf(:,:,:) = 0.0_wp … … 888 887 IF( cn_cfg == 'wad' ) THEN 889 888 ! 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 (:,:,:)889 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 890 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 891 ssh (:,:,Kmm) = ssh(:,:,Kbb) 892 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 893 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 894 ELSE 896 895 ! 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 896 ssh(:,:,Kmm) = -ssh_ref 897 ssh(:,:,Kbb) = -ssh_ref 898 899 DO_2D( 1, 1, 1, 1 ) 900 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 901 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 902 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 903 ENDIF 904 END_2D 910 905 ENDIF !If test case else 911 906 912 907 ! Adjust vertical metrics for all wad 913 908 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &909 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 910 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 911 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 912 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 913 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 914 915 DO_2D( 1, 1, 1, 1 ) 916 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 917 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 918 ENDIF 919 END_2D 927 920 ! 928 921 ELSE 929 922 ! 930 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n 931 CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) 932 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 923 ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 924 ! 925 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 926 ! 927 ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 933 928 ! 934 929 DO jk=1,jpk 935 e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &930 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 936 931 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 937 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t _b!= 0 on land points932 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(:,:,:,Kbb) != 0 on land points 938 933 END DO 939 e3t_n(:,:,:) = e3t_b(:,:,:) 940 sshn(:,:) = sshb(:,:) ! needed later for gde3w 941 !!$ e3t_n(:,:,:)=e3t_0(:,:,:) 942 !!$ e3t_b(:,:,:)=e3t_0(:,:,:) 934 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 935 ssh(:,:,Kmm) = ssh(:,:,Kbb) ! needed later for gde3w 943 936 ! 944 937 END IF ! end of ll_wd edits … … 958 951 ! ! all cases ! 959 952 ! ! --------- ! 960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t _b(:,:,:), ldxios = lwxios )961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t _n(:,:,:), ldxios = lwxios )953 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 954 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 962 955 ! ! ----------------------- ! 963 956 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 992 985 !!---------------------------------------------------------------------- 993 986 ! 994 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :995 987 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 996 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 997 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run 988 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' ) 998 989 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 999 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' , lwp)990 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) 1000 991 IF(lwm) WRITE ( numond, nam_vvl ) 1001 992 ! … … 1019 1010 WRITE(numout,*) ' rn_rst_e3t = 0.e0' 1020 1011 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1021 WRITE(numout,*) ' rn_lf_cutoff = 1.0/r dt'1012 WRITE(numout,*) ' rn_lf_cutoff = 1.0/rn_Dt' 1022 1013 ELSE 1023 1014 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t … … 1034 1025 ! 1035 1026 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1036 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1037 1027 ! 1038 1028 IF(lwp) THEN ! Print the choice … … 1050 1040 END SUBROUTINE dom_vvl_ctl 1051 1041 1042 #endif 1043 1052 1044 !!====================================================================== 1053 1045 END MODULE domvvl
Note: See TracChangeset
for help on using the changeset viewer.
