- Timestamp:
- 2013-04-09T18:34:38+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r3294 r3865 6 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! ----------------------------------------------------------------------9 #if defined key_vvl 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 !! vvl option includes z_star and z_tilde coordinates 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_vvl' variable volume 12 12 !!---------------------------------------------------------------------- 13 !! dom_vvl : defined coefficients to distribute ssh on each layers14 13 !!---------------------------------------------------------------------- 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 16 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 18 !! dom_vvl_rst : read/write restart file 19 !! dom_vvl_ctl : Check the vvl options 20 !!---------------------------------------------------------------------- 21 !! * Modules used 15 22 USE oce ! ocean dynamics and tracers 16 23 USE dom_oce ! ocean space and time domain 17 USE sbc_oce ! surface boundary condition: ocean 18 USE phycst ! physical constants 24 USE sbc_oce ! ocean surface boundary condition 19 25 USE in_out_manager ! I/O manager 26 USE iom ! I/O manager library 27 USE restart ! ocean restart 20 28 USE lib_mpp ! distributed memory computing library 21 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 26 34 PRIVATE 27 35 28 PUBLIC dom_vvl ! called by domain.F90 29 PUBLIC dom_vvl_2 ! called by domain.F90 30 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 31 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 36 !! * Routine accessibility 37 PUBLIC dom_vvl_init ! called by domain.F90 38 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_swp ! called by step.F90 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 42 !!* Namelist nam_vvl 43 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 44 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 45 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 46 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 47 ! ! conservation: not used yet 48 REAL(wp) :: rn_ahe3 = 0.e0 ! thickness diffusion coefficient 49 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 50 51 !! * Module variables 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 53 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 54 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors 55 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 56 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 33 57 34 58 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra … … 39 63 # include "vectopt_loop_substitute.h90" 40 64 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 4.0 , NEMO Consortium (2011)65 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 42 66 !! $Id$ 43 67 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 68 !!---------------------------------------------------------------------- 45 CONTAINS 69 70 CONTAINS 46 71 47 72 INTEGER FUNCTION dom_vvl_alloc() 48 73 !!---------------------------------------------------------------------- 49 !! *** ROUTINE dom_vvl_alloc *** 50 !!---------------------------------------------------------------------- 51 ! 52 ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) , & 53 & r2dt (jpk) , STAT=dom_vvl_alloc ) 54 ! 55 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 56 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 57 ! 74 !! *** FUNCTION dom_vvl_alloc *** 75 !!---------------------------------------------------------------------- 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 ALLOCATE( e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) , & 79 & un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , STAT = dom_vvl_alloc ) 80 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 81 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 82 ENDIF 83 IF( ln_vvl_ztilde ) THEN 84 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 85 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 86 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 87 ENDIF 88 58 89 END FUNCTION dom_vvl_alloc 59 90 60 91 61 SUBROUTINE dom_vvl 62 !!---------------------------------------------------------------------- 63 !! *** ROUTINE dom_vvl ***92 SUBROUTINE dom_vvl_init 93 !!---------------------------------------------------------------------- 94 !! *** ROUTINE dom_vvl_init *** 64 95 !! 65 !! ** Purpose : compute mu coefficients at t-, u-, v- and f-points to 66 !! spread ssh over the whole water column (scale factors) 67 !! set the before and now ssh at u- and v-points 68 !! (also f-point in now case) 69 !!---------------------------------------------------------------------- 70 ! 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 REAL(wp) :: zcoefu, zcoefv , zcoeff ! local scalars 73 REAL(wp) :: zvt , zvt_ip1, zvt_jp1, zvt_ip1jp1 ! - - 74 REAL(wp), POINTER, DIMENSION(:,:) :: zee_t, zee_u, zee_v, zee_f ! 2D workspace 75 !!---------------------------------------------------------------------- 76 ! 77 IF( nn_timing == 1 ) CALL timing_start('dom_vvl') 78 ! 79 CALL wrk_alloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 80 ! 81 IF(lwp) THEN 82 WRITE(numout,*) 83 WRITE(numout,*) 'dom_vvl : Variable volume initialization' 84 WRITE(numout,*) '~~~~~~~~ compute coef. used to spread ssh over each layers' 85 ENDIF 86 87 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' ) 88 89 fsdept(:,:,:) = gdept (:,:,:) 90 fsdepw(:,:,:) = gdepw (:,:,:) 91 fsde3w(:,:,:) = gdep3w(:,:,:) 92 fse3t (:,:,:) = e3t (:,:,:) 93 fse3u (:,:,:) = e3u (:,:,:) 94 fse3v (:,:,:) = e3v (:,:,:) 95 fse3f (:,:,:) = e3f (:,:,:) 96 fse3w (:,:,:) = e3w (:,:,:) 97 fse3uw(:,:,:) = e3uw (:,:,:) 98 fse3vw(:,:,:) = e3vw (:,:,:) 99 100 ! !== mu computation ==! 101 zee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 102 zee_u(:,:) = fse3u_0(:,:,1) 103 zee_v(:,:) = fse3v_0(:,:,1) 104 zee_f(:,:) = fse3f_0(:,:,1) 105 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 106 zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 107 zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 108 zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 109 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 110 zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 111 END DO 112 END DO 113 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 114 zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 115 zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 116 zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 117 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 118 zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 96 !! ** Purpose : Initialization of all scale factors, depths 97 !! and water column heights 98 !! 99 !! ** Method : - use restart file and/or initialize 100 !! - interpolate scale factors 101 !! 102 !! ** Action : - fse3t_(n/b) and e3t_t_(n/b) 103 !! - Regrid: fse3(u/v)_n 104 !! fse3(u/v)_b 105 !! fse3w_n 106 !! fse3(u/v)w_b 107 !! fse3(u/v)w_n 108 !! fsdept_n, fsdepw_n and fsde3w_n 109 !! - h(t/u/v)_0 110 !! - frq_rst_e3t and frq_rst_hdv 111 !! 112 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 113 !!---------------------------------------------------------------------- 114 USE phycst, ONLY : rpi 115 !! * Local declarations 116 INTEGER :: jk 117 !!---------------------------------------------------------------------- 118 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 119 120 IF(lwp) WRITE(numout,*) 121 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 122 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 123 124 ! choose vertical coordinate (z_star, z_tilde or layer) 125 ! ========================== 126 CALL dom_vvl_ctl 127 128 ! Allocate module arrays 129 ! ====================== 130 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 131 132 ! Read or initialize fse3t_(b/n), e3t_t_(b/n) and hdiv_lf (and e3t_a(jpk)) 133 ! ======================================================================== 134 CALL dom_vvl_rst( nit000, 'READ' ) 135 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 136 137 ! Reconstruction of all vertical scale factors at now and before time steps 138 ! ========================================================================= 139 ! Horizontal scale factor interpolations 140 ! -------------------------------------- 141 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 142 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 143 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 144 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 145 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 146 ! Vertical scale factor interpolations 147 ! ------------------------------------ 148 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 149 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 150 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 151 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 152 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 153 ! t- and w- points depth 154 ! ---------------------- 155 fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 156 fsdepw_n(:,:,1) = 0.e0 157 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 158 DO jk = 2, jpk 159 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 160 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 161 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 119 162 END DO 120 CALL lbc_lnk( zee_f, 'F', 1. ) ! lateral boundary condition on ee_f 121 ! 122 DO jk = 1, jpk ! mu coefficients 123 mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 124 muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk) ! U-point at T levels 125 muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 126 END DO 127 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 128 DO jj = 1, jpjm1 129 muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 130 END DO 131 muf(:,jpj,jk) = 0._wp 132 END DO 133 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition 134 135 136 hu_0(:,:) = 0.e0 ! Reference ocean depth at U- and V-points 163 ! Reference water column height at t-, u- and v- point 164 ! ---------------------------------------------------- 165 ht_0(:,:) = 0.e0 166 hu_0(:,:) = 0.e0 137 167 hv_0(:,:) = 0.e0 138 168 DO jk = 1, jpk 139 hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 140 hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 169 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 170 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 171 hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 141 172 END DO 142 143 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 144 DO ji = 1, jpim1 ! NO vector opt. 145 zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 146 zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 147 zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 173 174 ! Restoring frequencies for z_tilde coordinate 175 ! ============================================ 176 IF( ln_vvl_ztilde ) THEN 177 ! - ML - In the future, this should be tunable by the user (namelist) 178 frq_rst_e3t(:,:) = 2.e0_wp * rpi / ( 30.e0_wp * 86400.e0_wp ) 179 frq_rst_hdv(:,:) = 2.e0_wp * rpi / ( 5.e0_wp * 86400.e0_wp ) 180 ENDIF 181 182 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 183 184 END SUBROUTINE dom_vvl_init 185 186 187 SUBROUTINE dom_vvl_sf_nxt( kt ) 188 !!---------------------------------------------------------------------- 189 !! *** ROUTINE dom_vvl_sf_nxt *** 190 !! 191 !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, 192 !! tranxt and dynspg routines 193 !! 194 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 195 !! - z_tilde_case: after scale factor increment = 196 !! high frequency part of horizontal divergence 197 !! + retsoring towards the background grid 198 !! + thickness difusion 199 !! Then repartition of ssh INCREMENT proportionnaly 200 !! to the "baroclinic" level thickness. 201 !! 202 !! ** Action : - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case 203 !! - e3t_t_a: after increment of vertical scale factor 204 !! in z_tilde case 205 !! - fse3(t/u/v)_a 206 !! 207 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 208 !!---------------------------------------------------------------------- 209 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 210 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 211 !! * Arguments 212 INTEGER, INTENT( in ) :: kt ! time step 213 !! * Local declarations 214 INTEGER :: ji, jj, jk ! dummy loop indices 215 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 216 REAL(wp) :: z2dt ! temporary scalars 217 REAL(wp) :: z_def_max ! temporary scalar 218 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 219 !!---------------------------------------------------------------------- 220 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 221 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 222 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 223 224 IF(kt == nit000) THEN 225 IF(lwp) WRITE(numout,*) 226 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 227 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 228 ENDIF 229 230 ! ******************************* ! 231 ! After acale factors at t-points ! 232 ! ******************************* ! 233 234 ! ! ----------------- ! 235 IF( ln_vvl_zstar ) THEN ! z_star coordinate ! 236 ! ! ----------------- ! 237 238 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 239 DO jk = 1, jpkm1 240 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 241 END DO 242 243 ! ! --------------------------- ! 244 ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 245 ! ! --------------------------- ! 246 247 ! I - initialization 248 ! ================== 249 250 ! 1 - barotropic divergence 251 ! ------------------------- 252 zhdiv(:,:) = 0. 253 zht(:,:) = 0. 254 DO jk = 1, jpkm1 255 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 256 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 257 END DO 258 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 259 260 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 261 ! -------------------------------------------------- 262 IF( ln_vvl_ztilde ) THEN 263 IF( kt .GT. nit000 ) THEN 264 DO jk = 1, jpkm1 265 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 266 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 267 END DO 268 ENDIF 269 END IF 270 271 ! II - after z_tilde increments of vertical scale factors 272 ! ======================================================= 273 e3t_t_a(:,:,:) = 0.e0 ! e3t_t_a used to store tendency terms 274 275 ! 1 - High frequency divergence term 276 ! ---------------------------------- 277 IF( ln_vvl_ztilde ) THEN ! z_tilde case 278 DO jk = 1, jpkm1 279 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 280 END DO 281 ELSE ! layer case 282 DO jk = 1, jpkm1 283 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 284 END DO 285 END IF 286 287 ! 2 - Restoring term (z-tilde case only) 288 ! ------------------ 289 IF( ln_vvl_ztilde ) THEN 290 DO jk = 1, jpk 291 e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_rst_e3t(:,:) * e3t_t_b(:,:,jk) 292 END DO 293 END IF 294 295 ! 3 - Thickness diffusion term 296 ! ---------------------------- 297 zwu(:,:) = 0.e0 298 zwv(:,:) = 0.e0 299 ! a - first derivative: diffusive fluxes 300 DO jk = 1, jpkm1 301 DO jj = 1, jpjm1 302 DO ji = 1, fs_jpim1 ! vector opt. 303 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj ,jk) ) 304 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji ,jj+1,jk) ) 305 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 306 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 307 END DO 308 END DO 309 END DO 310 ! b - correction for last oceanic u-v points 311 DO jj = 1, jpj 312 DO ji = 1, jpi 313 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 314 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 315 END DO 316 END DO 317 ! c - second derivative: divergence of diffusive fluxes 318 DO jk = 1, jpkm1 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 322 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) ) & 323 & * r1_e12t(ji,jj) 324 END DO 325 END DO 326 END DO 327 ! d - thickness diffusion transport: boundary conditions 328 ! (stored for tracer advction and continuity equation) 329 CALL lbc_lnk( un_td , 'U' , -1.) 330 CALL lbc_lnk( vn_td , 'V' , -1.) 331 332 ! 4 - Time stepping of baroclinic scale factors 333 ! --------------------------------------------- 334 ! Leapfrog time stepping 335 ! ~~~~~~~~~~~~~~~~~~~~~~ 336 IF( neuler == 0 .AND. kt == nit000 ) THEN 337 z2dt = rdt 338 ELSE 339 z2dt = 2.e0 * rdt 340 ENDIF 341 CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. ) 342 e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:) 343 344 ! Maximum deformation control 345 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 346 ! - ML - Should perhaps be put in the namelist 347 z_def_max = 0.9e0 348 ze3t(:,:,jpk) = 0.e0 349 DO jk = 1, jpkm1 350 ze3t(:,:,jk) = e3t_t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 351 END DO 352 z_tmax = MAXVAL( ze3t(:,:,:) ) 353 z_tmin = MINVAL( ze3t(:,:,:) ) 354 ! - ML - test: for the moment, stop simulation for too large e3_t variations 355 IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN 356 ijk_max = MAXLOC( ze3t(:,:,:) ) 357 ijk_min = MINLOC( ze3t(:,:,:) ) 358 WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 359 WRITE(numout, *) 'at i, j, k=', ijk_max 360 WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 361 WRITE(numout, *) 'at i, j, k=', ijk_min 362 CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 363 ENDIF 364 ! - ML - end test 365 ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used 366 e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:), z_def_max * e3t_0(:,:,:) ) 367 e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), - z_def_max * e3t_0(:,:,:) ) 368 369 ! Add "tilda" part to the after scale factor 370 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 371 fse3t_a(:,:,:) = e3t_0(:,:,:) + e3t_t_a(:,:,:) 372 373 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 374 ! ================================================================================== 375 ! add e3t(n-1) "star" Asselin-filtered 376 DO jk = 1, jpkm1 377 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - e3t_t_b(:,:,jk) 378 END DO 379 ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n) 380 ! - ML - baroclinicity error should be better treated in the future 381 ! i.e. locally and not spread over the water column. 382 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 383 zht(:,:) = 0. 384 DO jk = 1, jpkm1 385 zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk) 386 END DO 387 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 388 DO jk = 1, jpkm1 389 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 390 END DO 391 392 ENDIF 393 394 IF( ln_vvl_dbg ) THEN ! - ML - test: control prints for debuging 395 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 396 WRITE(numout, *) 'kt =', kt 397 WRITE(numout, *) 'MAXVAL(abs(SUM(e3t_t_a))) =', & 398 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 399 END IF 400 zht(:,:) = 0.e0 401 DO jk = 1, jpkm1 402 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 403 END DO 404 WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', & 405 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 406 zht(:,:) = 0.e0 407 DO jk = 1, jpkm1 408 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 409 END DO 410 WRITE(numout, *) 'MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', & 411 & MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 412 END IF 413 414 ! *********************************** ! 415 ! After scale factors at u- v- points ! 416 ! *********************************** ! 417 418 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' ) 419 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' ) 420 421 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 422 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 423 424 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 425 426 END SUBROUTINE dom_vvl_sf_nxt 427 428 429 SUBROUTINE dom_vvl_sf_swp( kt ) 430 !!---------------------------------------------------------------------- 431 !! *** ROUTINE dom_vvl_sf_swp *** 432 !! 433 !! ** Purpose : compute time filter and swap of scale factors 434 !! compute all depths and related variables for next time step 435 !! write outputs and restart file 436 !! 437 !! ** Method : - swap of e3t with trick for volume/tracer conservation 438 !! - reconstruct scale factor at other grid points (interpolate) 439 !! - recompute depths and water height fields 440 !! 441 !! ** Action : - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step 442 !! - Recompute: 443 !! fse3(u/v)_b 444 !! fse3w_n 445 !! fse3(u/v)w_b 446 !! fse3(u/v)w_n 447 !! fsdept_n, fsdepw_n and fsde3w_n 448 !! h(u/v) and h(u/v)r 449 !! 450 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 451 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 452 !!---------------------------------------------------------------------- 453 !! * Arguments 454 INTEGER, INTENT( in ) :: kt ! time step 455 !! * Local declarations 456 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 457 INTEGER :: jk ! dummy loop indices 458 !!---------------------------------------------------------------------- 459 460 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 461 ! 462 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def ) 463 ! 464 IF( kt == nit000 ) THEN 465 IF(lwp) WRITE(numout,*) 466 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 467 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step' 468 ENDIF 469 ! 470 ! Time filter and swap of scale factors 471 ! ===================================== 472 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt. 473 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 474 IF( neuler == 0 .AND. kt == nit000 ) THEN 475 e3t_t_b(:,:,:) = e3t_t_n(:,:,:) 476 ELSE 477 e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * ( e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) ) 478 ENDIF 479 e3t_t_n(:,:,:) = e3t_t_a(:,:,:) 480 ENDIF 481 fse3t_n(:,:,:) = fse3t_a(:,:,:) 482 fse3u_n(:,:,:) = fse3u_a(:,:,:) 483 fse3v_n(:,:,:) = fse3v_a(:,:,:) 484 485 ! Compute all missing vertical scale factor and depths 486 ! ==================================================== 487 ! Horizontal scale factor interpolations 488 ! -------------------------------------- 489 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 490 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 491 ! Vertical scale factor interpolations 492 ! ------------------------------------ 493 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 494 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 495 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 496 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 497 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 498 ! t- and w- points depth 499 ! ---------------------- 500 fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1) 501 fsdepw_n(:,:,1) = 0.e0 502 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 503 DO jk = 2, jpk 504 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 505 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 506 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 507 END DO 508 ! Local depth and Inverse of the local depth of the water column at u- and v- points 509 ! ---------------------------------------------------------------------------------- 510 hu(:,:) = 0. 511 hv(:,:) = 0. 512 DO jk = 1, jpk 513 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 514 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 515 END DO 516 ! Inverse of the local depth 517 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) ) 518 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) ) 519 520 ! Write outputs 521 ! ============= 522 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 523 CALL iom_put( "e3t_n" , fse3t_n (:,:,:) ) 524 CALL iom_put( "dept_n" , fsde3w_n (:,:,:) ) 525 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 526 527 ! write restart file 528 ! ================== 529 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 530 ! 531 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 532 ! 533 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 534 535 END SUBROUTINE dom_vvl_sf_swp 536 537 538 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 539 !!--------------------------------------------------------------------- 540 !! *** ROUTINE dom_vvl__interpol *** 541 !! 542 !! ** Purpose : interpolate scale factors from one grid point to another 543 !! 544 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0) 545 !! - horizontal interpolation: grid cell surface averaging 546 !! - vertical interpolation: simple averaging 547 !!---------------------------------------------------------------------- 548 !! * Arguments 549 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 550 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 551 CHARACTER(LEN=1), INTENT( in ) :: pout ! grid point of out scale factors 552 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 553 !! * Local declarations 554 INTEGER :: ji, jj, jk ! dummy loop indices 555 !!---------------------------------------------------------------------- 556 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 557 SELECT CASE ( pout ) 558 ! ! ------------------------------------- ! 559 CASE( 'U' ) ! interpolation from T-point to U-point ! 560 ! ! ------------------------------------- ! 561 ! horizontal surface weighted interpolation 562 DO jk = 1, jpk 563 DO jj = 1, jpjm1 564 DO ji = 1, fs_jpim1 ! vector opt. 565 pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * r1_e12u(ji,jj) & 566 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 567 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 568 END DO 569 END DO 570 END DO 571 ! boundary conditions 572 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. ) 573 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 574 ! ! ------------------------------------- ! 575 CASE( 'V' ) ! interpolation from T-point to V-point ! 576 ! ! ------------------------------------- ! 577 ! horizontal surface weighted interpolation 578 DO jk = 1, jpk 579 DO jj = 1, jpjm1 580 DO ji = 1, fs_jpim1 ! vector opt. 581 pe3_out(ji,jj,jk) = 0.5 * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 582 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 583 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 584 END DO 585 END DO 586 END DO 587 ! boundary conditions 588 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. ) 589 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 590 ! ! ------------------------------------- ! 591 CASE( 'F' ) ! interpolation from U-point to F-point ! 592 ! ! ------------------------------------- ! 593 ! horizontal surface weighted interpolation 594 DO jk = 1, jpk 595 DO jj = 1, jpjm1 596 DO ji = 1, fs_jpim1 ! vector opt. 597 pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 598 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 599 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 600 END DO 601 END DO 602 END DO 603 ! boundary conditions 604 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. ) 605 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 606 ! ! ------------------------------------- ! 607 CASE( 'W' ) ! interpolation from T-point to W-point ! 608 ! ! ------------------------------------- ! 609 ! vertical simple interpolation 610 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 611 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 612 DO jk = 2, jpk 613 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 614 & + 0.5 * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 615 END DO 616 ! ! -------------------------------------- ! 617 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 618 ! ! -------------------------------------- ! 619 ! vertical simple interpolation 620 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 621 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 622 DO jk = 2, jpk 623 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 624 & + 0.5 * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 625 END DO 626 ! ! -------------------------------------- ! 627 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 628 ! ! -------------------------------------- ! 629 ! vertical simple interpolation 630 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 631 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 632 DO jk = 2, jpk 633 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 634 & + 0.5 * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 635 END DO 636 END SELECT 637 638 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 639 640 END SUBROUTINE dom_vvl_interpol 641 642 643 SUBROUTINE dom_vvl_rst( kt, cdrw ) 644 !!--------------------------------------------------------------------- 645 !! *** ROUTINE dom_vvl_rst *** 646 !! 647 !! ** Purpose : Read or write VVL file in restart file 648 !! 649 !! ** Method : use of IOM library 650 !! if the restart does not contain vertical scale factors, 651 !! they are set to the _0 values 652 !! if the restart does not contain vertical scale factors increments (z_tilde), 653 !! they are set to 0. 654 !!---------------------------------------------------------------------- 655 !! * Arguments 656 INTEGER , INTENT(in) :: kt ! ocean time-step 657 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 658 !! * Local declarations 659 INTEGER :: id1, id2, id3, id4, id5 ! local integers 660 !!---------------------------------------------------------------------- 661 ! 662 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst') 663 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 664 ! ! =============== 665 IF( ln_rstart ) THEN !* Read the restart file 666 id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 667 id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 668 id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. ) 669 id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. ) 670 id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. ) 671 ! ! --------- ! 672 ! ! all cases ! 673 ! ! --------- ! 674 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 675 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 676 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 677 IF( neuler == 0 ) THEN 678 fse3t_b(:,:,:) = fse3t_n(:,:,:) 679 ENDIF 680 ELSE ! one at least array is missing 681 CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' ) 682 ENDIF 683 ! ! ----------- ! 684 IF( ln_vvl_zstar ) THEN ! z_star case ! 685 ! ! ----------- ! 686 IF( MIN( id3, id4 ) > 0 ) THEN 687 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 688 ENDIF 689 ! ! ----------------------- ! 690 ELSE ! z_tilde and layer cases ! 691 ! ! ----------------------- ! 692 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 693 CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) ) 694 CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) ) 695 ELSE ! one at least array is missing 696 e3t_t_b(:,:,:) = 0.e0 697 e3t_t_n(:,:,:) = 0.e0 698 ENDIF 699 ! ! ------------ ! 700 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 701 ! ! ------------ ! 702 IF( id5 > 0 ) THEN ! required array exists 703 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 704 ELSE ! array is missing 705 hdiv_lf(:,:,:) = 0.e0 706 ENDIF 707 ENDIF 708 ENDIF 148 709 ! 149 zvt = e1e2t(ji ,jj ) * sshb(ji ,jj ) ! before fields 150 zvt_ip1 = e1e2t(ji+1,jj ) * sshb(ji+1,jj ) 151 zvt_jp1 = e1e2t(ji ,jj+1) * sshb(ji ,jj+1) 152 sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 153 sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 154 ! 155 zvt = e1e2t(ji ,jj ) * sshn(ji ,jj ) ! now fields 156 zvt_ip1 = e1e2t(ji+1,jj ) * sshn(ji+1,jj ) 157 zvt_jp1 = e1e2t(ji ,jj+1) * sshn(ji ,jj+1) 158 zvt_ip1jp1 = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 159 sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 160 sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 161 sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 162 END DO 163 END DO 164 CALL lbc_lnk( sshu_n, 'U', 1. ) ; CALL lbc_lnk( sshu_b, 'U', 1. ) ! lateral boundary conditions 165 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 166 CALL lbc_lnk( sshf_n, 'F', 1. ) 167 ! 168 CALL wrk_dealloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 169 ! 170 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl') 171 ! 172 END SUBROUTINE dom_vvl 173 174 175 SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 176 !!---------------------------------------------------------------------- 177 !! *** ROUTINE dom_vvl_2 *** 178 !! 179 !! ** Purpose : compute the vertical scale factors at u- and v-points 180 !! in variable volume case. 181 !! 182 !! ** Method : In variable volume case (non linear sea surface) the 183 !! the vertical scale factor at velocity points is computed 184 !! as the average of the cell surface weighted e3t. 185 !! It uses the sea surface heigth so it have to be initialized 186 !! after ssh is read/set 187 !!---------------------------------------------------------------------- 188 INTEGER , INTENT(in ) :: kt ! ocean time-step index 189 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3u_b, pe3v_b ! before vertical scale factor at u- & v-pts 190 ! 191 INTEGER :: ji, jj, jk ! dummy loop indices 192 INTEGER :: iku, ikv ! local integers 193 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 194 REAL(wp) :: zvt ! local scalars 195 !!---------------------------------------------------------------------- 196 ! 197 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_2') 198 ! 199 IF( lwp .AND. kt == nit000 ) THEN 710 ELSE !* Initialize at "rest" 711 fse3t_b(:,:,:) = e3t_0(:,:,:) 712 fse3t_n(:,:,:) = e3t_0(:,:,:) 713 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 714 e3t_t_b(:,:,:) = 0.e0 715 e3t_t_n(:,:,:) = 0.e0 716 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0 717 END IF 718 ENDIF 719 720 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 721 ! ! =================== 722 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 723 ! ! --------- ! 724 ! ! all cases ! 725 ! ! --------- ! 726 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 727 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) 728 ! ! ----------------------- ! 729 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 730 ! ! ----------------------- ! 731 CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) ) 732 CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) ) 733 END IF 734 ! ! -------------! 735 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 736 ! ! ------------ ! 737 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 738 ENDIF 739 740 ENDIF 741 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 742 743 END SUBROUTINE dom_vvl_rst 744 745 746 SUBROUTINE dom_vvl_ctl 747 !!--------------------------------------------------------------------- 748 !! *** ROUTINE dom_vvl_ctl *** 749 !! 750 !! ** Purpose : Control the consistency between namelist options 751 !! for vertical coordinate 752 !!---------------------------------------------------------------------- 753 INTEGER :: ioptio 754 755 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, rn_ahe3, ln_vvl_dbg! , ln_vvl_kepe 756 !!---------------------------------------------------------------------- 757 758 REWIND ( numnam ) ! Read Namelist nam_vvl : vertical coordinate 759 READ ( numnam, nam_vvl ) 760 761 IF(lwp) THEN ! Namelist print 200 762 WRITE(numout,*) 201 WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 202 WRITE(numout,*) '~~~~~~~~~ ' 203 pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 204 pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 205 ENDIF 206 207 DO jk = 1, jpkm1 ! set the before scale factors at u- & v-points 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 210 zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 211 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 212 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 213 END DO 214 END DO 215 END DO 216 217 ! Correct scale factors at locations that have been individually modified in domhgr 218 ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute 219 ! scale factors ignoring the modified metric. 220 ! ! ===================== 221 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 222 ! ! ===================== 223 IF( nn_cla == 0 ) THEN 224 ! 225 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified) 226 ij0 = 102 ; ij1 = 102 227 DO jk = 1, jpkm1 ! set the before scale factors at u-points 228 DO jj = mj0(ij0), mj1(ij1) 229 DO ji = mi0(ii0), mi1(ii1) 230 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 231 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 232 END DO 233 END DO 234 END DO 235 ! 236 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified) 237 ij0 = 88 ; ij1 = 88 238 DO jk = 1, jpkm1 ! set the before scale factors at u-points 239 DO jj = mj0(ij0), mj1(ij1) 240 DO ji = mi0(ii0), mi1(ii1) 241 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 242 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 243 END DO 244 END DO 245 END DO 246 DO jk = 1, jpkm1 ! set the before scale factors at v-points 247 DO jj = mj0(ij0), mj1(ij1) 248 DO ji = mi0(ii0), mi1(ii1) 249 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 250 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 251 END DO 252 END DO 253 END DO 254 ENDIF 255 256 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified) 257 ij0 = 116 ; ij1 = 116 258 DO jk = 1, jpkm1 ! set the before scale factors at u-points 259 DO jj = mj0(ij0), mj1(ij1) 260 DO ji = mi0(ii0), mi1(ii1) 261 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 262 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 263 END DO 264 END DO 265 END DO 266 ! 267 ENDIF 268 ! ! ===================== 269 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 270 ! ! ===================== 271 272 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified) 273 ij0 = 200 ; ij1 = 200 274 DO jk = 1, jpkm1 ! set the before scale factors at u-points 275 DO jj = mj0(ij0), mj1(ij1) 276 DO ji = mi0(ii0), mi1(ii1) 277 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 278 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 279 END DO 280 END DO 281 END DO 282 283 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified) 284 ij0 = 208 ; ij1 = 208 285 DO jk = 1, jpkm1 ! set the before scale factors at u-points 286 DO jj = mj0(ij0), mj1(ij1) 287 DO ji = mi0(ii0), mi1(ii1) 288 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 289 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 290 END DO 291 END DO 292 END DO 293 294 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified) 295 ij0 = 124 ; ij1 = 125 296 DO jk = 1, jpkm1 ! set the before scale factors at v-points 297 DO jj = mj0(ij0), mj1(ij1) 298 DO ji = mi0(ii0), mi1(ii1) 299 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 300 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 301 END DO 302 END DO 303 END DO 304 305 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 306 ij0 = 124 ; ij1 = 125 307 DO jk = 1, jpkm1 ! set the before scale factors at v-points 308 DO jj = mj0(ij0), mj1(ij1) 309 DO ji = mi0(ii0), mi1(ii1) 310 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 311 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 312 END DO 313 END DO 314 END DO 315 316 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified) 317 ij0 = 124 ; ij1 = 125 318 DO jk = 1, jpkm1 ! set the before scale factors at v-points 319 DO jj = mj0(ij0), mj1(ij1) 320 DO ji = mi0(ii0), mi1(ii1) 321 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 322 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 323 END DO 324 END DO 325 END DO 326 327 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified) 328 ij0 = 124 ; ij1 = 125 329 DO jk = 1, jpkm1 ! set the before scale factors at v-points 330 DO jj = mj0(ij0), mj1(ij1) 331 DO ji = mi0(ii0), mi1(ii1) 332 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 333 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 334 END DO 335 END DO 336 END DO 337 338 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified) 339 ij0 = 141 ; ij1 = 142 340 DO jk = 1, jpkm1 ! set the before scale factors at v-points 341 DO jj = mj0(ij0), mj1(ij1) 342 DO ji = mi0(ii0), mi1(ii1) 343 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 344 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 345 END DO 346 END DO 347 END DO 348 349 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified) 350 ij0 = 141 ; ij1 = 142 351 DO jk = 1, jpkm1 ! set the before scale factors at v-points 352 DO jj = mj0(ij0), mj1(ij1) 353 DO ji = mi0(ii0), mi1(ii1) 354 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 355 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 356 END DO 357 END DO 358 END DO 359 360 ! 361 ENDIF 362 ! ! ====================== 363 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 364 ! ! ====================== 365 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified) 366 ij0 = 327 ; ij1 = 327 367 DO jk = 1, jpkm1 ! set the before scale factors at u-points 368 DO jj = mj0(ij0), mj1(ij1) 369 DO ji = mi0(ii0), mi1(ii1) 370 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 371 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 372 END DO 373 END DO 374 END DO 375 ! 376 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u was modified) 377 ij0 = 343 ; ij1 = 343 378 DO jk = 1, jpkm1 ! set the before scale factors at u-points 379 DO jj = mj0(ij0), mj1(ij1) 380 DO ji = mi0(ii0), mi1(ii1) 381 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 382 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 383 END DO 384 END DO 385 END DO 386 ! 387 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified) 388 ij0 = 232 ; ij1 = 232 389 DO jk = 1, jpkm1 ! set the before scale factors at u-points 390 DO jj = mj0(ij0), mj1(ij1) 391 DO ji = mi0(ii0), mi1(ii1) 392 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 393 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 394 END DO 395 END DO 396 END DO 397 ! 398 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified) 399 ij0 = 232 ; ij1 = 232 400 DO jk = 1, jpkm1 ! set the before scale factors at u-points 401 DO jj = mj0(ij0), mj1(ij1) 402 DO ji = mi0(ii0), mi1(ii1) 403 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 404 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 405 END DO 406 END DO 407 END DO 408 ! 409 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified) 410 ij0 = 270 ; ij1 = 270 411 DO jk = 1, jpkm1 ! set the before scale factors at u-points 412 DO jj = mj0(ij0), mj1(ij1) 413 DO ji = mi0(ii0), mi1(ii1) 414 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 415 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 416 END DO 417 END DO 418 END DO 419 ! 420 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified) 421 ij0 = 232 ; ij1 = 233 422 DO jk = 1, jpkm1 ! set the before scale factors at v-points 423 DO jj = mj0(ij0), mj1(ij1) 424 DO ji = mi0(ii0), mi1(ii1) 425 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 426 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 427 END DO 428 END DO 429 END DO 430 ! 431 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified) 432 ij0 = 276 ; ij1 = 276 433 DO jk = 1, jpkm1 ! set the before scale factors at v-points 434 DO jj = mj0(ij0), mj1(ij1) 435 DO ji = mi0(ii0), mi1(ii1) 436 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 437 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 438 END DO 439 END DO 440 END DO 441 ! 442 ENDIF 443 ! End of individual corrections to scale factors 444 445 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 446 DO jj = 2, jpjm1 447 DO ji = fs_2, fs_jpim1 448 iku = mbku(ji,jj) 449 ikv = mbkv(ji,jj) 450 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 451 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 452 END DO 453 END DO 454 ENDIF 455 456 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos 457 pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 458 CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 459 CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 460 pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:) ! recover the full scale factor 461 pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 462 ! 463 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_2') 464 ! 465 END SUBROUTINE dom_vvl_2 466 467 #else 468 !!---------------------------------------------------------------------- 469 !! Default option : Empty routine 470 !!---------------------------------------------------------------------- 471 CONTAINS 472 SUBROUTINE dom_vvl 473 END SUBROUTINE dom_vvl 474 SUBROUTINE dom_vvl_2(kdum, pudum, pvdum ) 475 USE par_kind 476 INTEGER , INTENT(in ) :: kdum 477 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pudum, pvdum 478 END SUBROUTINE dom_vvl_2 479 #endif 763 WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 764 WRITE(numout,*) '~~~~~~~~~~~' 765 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 766 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar 767 WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde 768 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 769 ! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation' 770 ! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe 771 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient' 772 WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3 773 WRITE(numout,*) ' Namelist nam_vvl : debug prints' 774 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 775 ENDIF 776 777 ioptio = 0 ! Parameter control 778 IF( ln_vvl_zstar ) ioptio = ioptio + 1 779 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 780 IF( ln_vvl_layer ) ioptio = ioptio + 1 781 782 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 783 784 IF(lwp) THEN ! Print the choice 785 WRITE(numout,*) 786 IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used' 787 IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used' 788 IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used' 789 ! - ML - Option not developed yet 790 ! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used' 791 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 792 ENDIF 793 794 END SUBROUTINE dom_vvl_ctl 480 795 481 796 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
