- Timestamp:
- 2015-10-31T08:40:45+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5836 r5845 20 20 !!---------------------------------------------------------------------- 21 21 USE oce ! ocean dynamics and tracers 22 USE phycst ! physical constant 22 23 USE dom_oce ! ocean space and time domain 23 24 USE sbc_oce ! ocean surface boundary condition 25 USE restart ! ocean restart 26 ! 24 27 USE in_out_manager ! I/O manager 25 28 USE iom ! I/O manager library 26 USE restart ! ocean restart27 29 USE lib_mpp ! distributed memory computing library 28 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 60 62 61 63 !! * Substitutions 62 # include "domzgr_substitute.h90"63 64 # include "vectopt_loop_substitute.h90" 64 65 !!---------------------------------------------------------------------- … … 67 68 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 68 69 !!---------------------------------------------------------------------- 69 70 70 CONTAINS 71 71 … … 81 81 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 82 82 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 83 un_td = 0. 0_wp84 vn_td = 0. 0_wp83 un_td = 0._wp 84 vn_td = 0._wp 85 85 ENDIF 86 86 IF( ln_vvl_ztilde ) THEN … … 89 89 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 90 90 ENDIF 91 91 ! 92 92 END FUNCTION dom_vvl_alloc 93 93 … … 103 103 !! - interpolate scale factors 104 104 !! 105 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)106 !! - Regrid: fse3(u/v)_n107 !! fse3(u/v)_b108 !! fse3w_n109 !! fse3(u/v)w_b110 !! fse3(u/v)w_n111 !! fsdept_n, fsdepw_n and fsde3w_n105 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 106 !! - Regrid: e3(u/v)_n 107 !! e3(u/v)_b 108 !! e3w_n 109 !! e3(u/v)w_b 110 !! e3(u/v)w_n 111 !! gdept_n, gdepw_n and gde3w_n 112 112 !! - h(t/u/v)_0 113 113 !! - frq_rst_e3t and frq_rst_hdv … … 115 115 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 116 !!---------------------------------------------------------------------- 117 USE phycst, ONLY : rpi, rsmall, rad 118 !! * Local declarations 119 INTEGER :: ji,jj,jk 117 INTEGER :: ji, jj, jk 120 118 INTEGER :: ii0, ii1, ij0, ij1 121 119 REAL(wp):: zcoef 122 120 !!---------------------------------------------------------------------- 121 ! 123 122 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 123 ! 125 124 IF(lwp) WRITE(numout,*) 126 125 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 127 126 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 128 129 ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! ========================== 131 CALL dom_vvl_ctl 132 133 ! Allocate module arrays 134 ! ====================== 127 ! 128 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 129 ! 130 ! ! Allocate module arrays 135 131 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 136 137 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 138 ! ============================================================================ 132 ! 133 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 139 134 CALL dom_vvl_rst( nit000, 'READ' ) 140 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 141 142 ! Reconstruction of all vertical scale factors at now and before time steps 143 ! ============================================================================= 144 ! Horizontal scale factor interpolations 145 ! -------------------------------------- 146 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 147 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 148 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 149 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 150 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 151 ! Vertical scale factor interpolations 152 ! ------------------------------------ 153 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 154 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 155 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 156 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 157 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 158 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 159 ! t- and w- points depth 160 ! ---------------------- 161 ! set the isf depth as it is in the initial step 162 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 163 fsdepw_n(:,:,1) = 0.0_wp 164 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 165 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 166 fsdepw_b(:,:,1) = 0.0_wp 167 168 DO jk = 2, jpk 135 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 136 ! 137 ! !== Set of all other vertical scale factors ==! (now and before) 138 ! ! Horizontal interpolation of e3t 139 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 140 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 141 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 142 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 143 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 144 ! ! Vertical interpolation of e3t,u,v 145 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 146 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 147 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 148 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 149 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 150 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 151 ! 152 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 153 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 154 gdepw_n(:,:,1) = 0.0_wp 155 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 156 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 157 gdepw_b(:,:,1) = 0.0_wp 158 DO jk = 2, jpk ! vertical sum 169 159 DO jj = 1,jpj 170 160 DO ji = 1,jpi 171 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 172 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 173 ! 0.5 where jk = mikt 174 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 175 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 176 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 177 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 178 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 179 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 180 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 181 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 161 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 162 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 163 ! ! 0.5 where jk = mikt 164 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 165 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 166 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 167 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 168 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 169 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 170 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 171 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 172 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 182 173 END DO 183 174 END DO 184 175 END DO 185 186 ! Before depth and Inverse of the local depth of the water column at u- and v- points 187 ! ----------------------------------------------------------------------------------- 188 hu_b(:,:) = 0. 189 hv_b(:,:) = 0. 190 DO jk = 1, jpkm1 191 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 192 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 176 ! 177 ! !== thickness of the water column !! (ocean portion only) 178 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 179 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 180 hv_b(:,:) = e3u_b(:,:,1) * vmask(:,:,1) 181 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 182 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 183 DO jk = 2, jpkm1 184 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 185 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 186 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 187 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 188 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 193 189 END DO 194 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 195 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 196 197 ! Restoring frequencies for z_tilde coordinate 198 ! ============================================ 190 ! 191 ! !== inverse of water column thickness ==! (u- and v- points) 192 r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) ! _i mask due to ISF 193 r1_hu_n(:,:) = umask_i(:,:) / ( hu_n(:,:) + 1._wp - umask_i(:,:) ) 194 r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 195 r1_hv_n(:,:) = vmask_i(:,:) / ( hv_n(:,:) + 1._wp - vmask_i(:,:) ) 196 197 ! !== z_tilde coordinate case ==! (Restoring frequencies) 199 198 IF( ln_vvl_ztilde ) THEN 200 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 201 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 202 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 203 IF( ln_vvl_ztilde_as_zstar ) THEN 204 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 205 frq_rst_e3t(:,:) = 0.0_wp 206 frq_rst_hdv(:,:) = 1.0_wp / rdt 207 ENDIF 208 IF ( ln_vvl_zstar_at_eqtor ) THEN 199 !!gm : idea: add here a READ in a file of custumized restoring frequency 200 ! ! Values in days provided via the namelist 201 ! ! use rsmall to avoid possible division by zero errors with faulty settings 202 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 203 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 204 ! 205 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 206 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 207 frq_rst_hdv(:,:) = 1._wp / rdt 208 ENDIF 209 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 209 210 DO jj = 1, jpj 210 211 DO ji = 1, jpi 212 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 211 213 IF( ABS(gphit(ji,jj)) >= 6.) THEN 212 214 ! values outside the equatorial band and transition zone (ztilde) 213 215 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 214 216 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 215 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 217 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 216 218 ! values inside the equatorial band (ztilde as zstar) 217 219 frq_rst_e3t(ji,jj) = 0.0_wp 218 220 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 219 ELSE 220 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)221 ELSE ! transition band (2.5 to 6 degrees N/S) 222 ! ! (linearly transition from z-tilde to z-star) 221 223 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 222 224 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 229 231 END DO 230 232 END DO 231 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 232 ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2233 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 234 ii0 = 103 ; ii1 = 111 233 235 ij0 = 128 ; ij1 = 135 ; 234 236 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 237 239 ENDIF 238 240 ENDIF 239 241 ! 240 242 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 241 243 ! 242 244 END SUBROUTINE dom_vvl_init 243 245 … … 261 263 !! - tilde_e3t_a: after increment of vertical scale factor 262 264 !! in z_tilde case 263 !! - fse3(t/u/v)_a265 !! - e3(t/u/v)_a 264 266 !! 265 267 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. … … 277 279 LOGICAL :: ll_do_bclinic ! temporary logical 278 280 !!---------------------------------------------------------------------- 281 ! 279 282 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 280 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 281 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 283 ! 284 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv ) 285 CALL wrk_alloc( jpi,jpj,jpk, ze3t ) 282 286 283 287 IF(kt == nit000) THEN … … 289 293 ll_do_bclinic = .TRUE. 290 294 IF( PRESENT(kcall) ) THEN 291 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.295 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 292 296 ENDIF 293 297 … … 295 299 ! After acale factors at t-points ! 296 300 ! ******************************* ! 297 298 301 ! ! --------------------------------------------- ! 299 302 ! ! z_star coordinate and barotropic z-tilde part ! 300 303 ! ! --------------------------------------------- ! 301 304 ! 302 305 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 303 306 DO jk = 1, jpkm1 304 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)305 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)307 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 308 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 306 309 END DO 307 310 ! 308 311 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 309 312 ! ! ------baroclinic part------ ! … … 314 317 ! 1 - barotropic divergence 315 318 ! ------------------------- 316 zhdiv(:,:) = 0. 317 zht(:,:) = 0. 319 zhdiv(:,:) = 0._wp 320 zht(:,:) = 0._wp 318 321 DO jk = 1, jpkm1 319 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)320 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)322 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 323 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 321 324 END DO 322 325 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 325 328 ! -------------------------------------------------- 326 329 IF( ln_vvl_ztilde ) THEN 327 IF( kt .GT.nit000 ) THEN330 IF( kt > nit000 ) THEN 328 331 DO jk = 1, jpkm1 329 332 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 330 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )333 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 331 334 END DO 332 335 ENDIF 333 END 336 ENDIF 334 337 335 338 ! II - after z_tilde increments of vertical scale factors 336 339 ! ======================================================= 337 tilde_e3t_a(:,:,:) = 0. 0_wp ! tilde_e3t_a used to store tendency terms340 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 338 341 339 342 ! 1 - High frequency divergence term … … 341 344 IF( ln_vvl_ztilde ) THEN ! z_tilde case 342 345 DO jk = 1, jpkm1 343 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )346 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 344 347 END DO 345 348 ELSE ! layer case 346 349 DO jk = 1, jpkm1 347 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)348 END DO 349 END 350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 351 END DO 352 ENDIF 350 353 351 354 ! 2 - Restoring term (z-tilde case only) … … 355 358 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 356 359 END DO 357 END 360 ENDIF 358 361 359 362 ! 3 - Thickness diffusion term 360 363 ! ---------------------------- 361 zwu(:,:) = 0.0_wp 362 zwv(:,:) = 0.0_wp 363 ! a - first derivative: diffusive fluxes 364 DO jk = 1, jpkm1 364 zwu(:,:) = 0._wp 365 zwv(:,:) = 0._wp 366 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 365 367 DO jj = 1, jpjm1 366 368 DO ji = 1, fs_jpim1 ! vector opt. … … 374 376 END DO 375 377 END DO 376 ! b - correction for last oceanic u-v points 377 DO jj = 1, jpj 378 DO jj = 1, jpj ! b - correction for last oceanic u-v points 378 379 DO ji = 1, jpi 379 380 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 381 382 END DO 382 383 END DO 383 ! c - second derivative: divergence of diffusive fluxes 384 DO jk = 1, jpkm1 384 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 385 385 DO jj = 2, jpjm1 386 386 DO ji = fs_2, fs_jpim1 ! vector opt. … … 391 391 END DO 392 392 END DO 393 ! d - thickness diffusion transport: boundary conditions394 ! (stored for tracer advction and continuity equation)393 ! ! d - thickness diffusion transport: boundary conditions 394 ! (stored for tracer advction and continuity equation) 395 395 CALL lbc_lnk( un_td , 'U' , -1._wp) 396 396 CALL lbc_lnk( vn_td , 'V' , -1._wp) … … 410 410 ! Maximum deformation control 411 411 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 412 ze3t(:,:,jpk) = 0. 0_wp412 ze3t(:,:,jpk) = 0._wp 413 413 DO jk = 1, jpkm1 414 414 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 462 462 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 463 463 DO jk = 1, jpkm1 464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 465 465 END DO 466 466 … … 470 470 ! ! ---baroclinic part--------- ! 471 471 DO jk = 1, jpkm1 472 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)472 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 473 473 END DO 474 474 ENDIF … … 485 485 zht(:,:) = 0.0_wp 486 486 DO jk = 1, jpkm1 487 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)487 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 488 488 END DO 489 489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 490 490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM( fse3t_n))) =', z_tmax491 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 492 492 ! 493 493 zht(:,:) = 0.0_wp 494 494 DO jk = 1, jpkm1 495 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)495 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 496 496 END DO 497 497 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 498 498 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 499 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM( fse3t_a))) =', z_tmax499 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 500 500 ! 501 501 zht(:,:) = 0.0_wp 502 502 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)503 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 504 504 END DO 505 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 506 506 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM( fse3t_b))) =', z_tmax507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 508 508 ! 509 509 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 524 524 ! *********************************** ! 525 525 526 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )527 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )526 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 527 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 528 528 529 529 ! *********************************** ! … … 531 531 ! *********************************** ! 532 532 533 hu_a(:,:) = 0._wp ! Ocean depth at U-points534 hv_a(:,:) = 0._wp ! Ocean depth at V-points535 DO jk = 1, jpkm1536 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)537 hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)533 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 534 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 535 DO jk = 2, jpkm1 536 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 537 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 538 538 END DO 539 539 ! ! Inverse of the local depth 540 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 541 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 542 543 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 544 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 545 540 !!gm BUG ? don't understand the use of umask_i here ..... 541 r1_hu_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 542 r1_hv_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 543 ! 544 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) 545 CALL wrk_dealloc( jpi,jpj,jpk, ze3t ) 546 ! 546 547 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 547 548 ! 548 549 END SUBROUTINE dom_vvl_sf_nxt 549 550 … … 561 562 !! - recompute depths and water height fields 562 563 !! 563 !! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step564 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 564 565 !! - Recompute: 565 !! fse3(u/v)_b566 !! fse3w_n567 !! fse3(u/v)w_b568 !! fse3(u/v)w_n569 !! fsdept_n, fsdepw_n and fsde3w_n566 !! e3(u/v)_b 567 !! e3w_n 568 !! e3(u/v)w_b 569 !! e3(u/v)w_n 570 !! gdept_n, gdepw_n and gde3w_n 570 571 !! h(u/v) and h(u/v)r 571 572 !! … … 573 574 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 574 575 !!---------------------------------------------------------------------- 575 !! * Arguments 576 INTEGER, INTENT( in ) :: kt ! time step 577 !! * Local declarations 578 INTEGER :: ji,jj,jk ! dummy loop indices 579 REAL(wp) :: zcoef 576 INTEGER, INTENT( in ) :: kt ! time step 577 ! 578 INTEGER :: ji, jj, jk ! dummy loop indices 579 REAL(wp) :: zcoef ! local scalar 580 580 !!---------------------------------------------------------------------- 581 581 … … 590 590 ! Time filter and swap of scale factors 591 591 ! ===================================== 592 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.592 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 593 593 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 594 594 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 600 600 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 601 601 ENDIF 602 fsdept_b(:,:,:) = fsdept_n(:,:,:)603 fsdepw_b(:,:,:) = fsdepw_n(:,:,:)604 605 fse3t_n(:,:,:) = fse3t_a(:,:,:)606 fse3u_n(:,:,:) = fse3u_a(:,:,:)607 fse3v_n(:,:,:) = fse3v_a(:,:,:)602 gdept_b(:,:,:) = gdept_n(:,:,:) 603 gdepw_b(:,:,:) = gdepw_n(:,:,:) 604 605 e3t_n(:,:,:) = e3t_a(:,:,:) 606 e3u_n(:,:,:) = e3u_a(:,:,:) 607 e3v_n(:,:,:) = e3v_a(:,:,:) 608 608 609 609 ! Compute all missing vertical scale factor and depths … … 611 611 ! Horizontal scale factor interpolations 612 612 ! -------------------------------------- 613 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt613 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 614 614 ! - JC - hu_b, hv_b, hur_b, hvr_b also 615 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 615 616 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 617 616 618 ! Vertical scale factor interpolations 617 619 ! ------------------------------------ 618 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )619 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )620 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )621 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' )622 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )623 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )620 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 621 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 622 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 623 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 624 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 625 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 624 626 ! t- and w- points depth 625 627 ! ---------------------- 626 628 ! set the isf depth as it is in the initial step 627 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)628 fsdepw_n(:,:,1) = 0.0_wp629 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)629 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 630 gdepw_n(:,:,1) = 0.0_wp 631 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 630 632 631 633 DO jk = 2, jpk … … 635 637 ! 1 for jk = mikt 636 638 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 637 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)638 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &639 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))640 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)639 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 640 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 641 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 642 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 641 643 END DO 642 644 END DO 643 645 END DO 644 646 645 ! Local depth and Inverse of the local depth of the water column at u- and v- points646 647 ! ---------------------------------------------------------------------------------- 647 hu (:,:) = hu_a (:,:) 648 hv (:,:) = hv_a (:,:) 649 650 ! Inverse of the local depth 651 hur(:,:) = hur_a(:,:) 652 hvr(:,:) = hvr_a(:,:) 653 654 ! Local depth of the water column at t- points 655 ! -------------------------------------------- 656 ht(:,:) = 0. 657 DO jk = 1, jpkm1 658 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 648 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 649 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 650 ! 651 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 652 DO jk = 2, jpkm1 653 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 659 654 END DO 660 655 661 656 ! Write outputs 662 657 ! ============= 663 CALL iom_put( "e3t" , fse3t_n (:,:,:) )664 CALL iom_put( "e3u" , fse3u_n (:,:,:) )665 CALL iom_put( "e3v" , fse3v_n (:,:,:) )666 CALL iom_put( "e3w" , fse3w_n (:,:,:) )667 CALL iom_put( "tpt_dep" , fsde3w_n(:,:,:) )658 CALL iom_put( "e3t" , e3t_n (:,:,:) ) 659 CALL iom_put( "e3u" , e3u_n (:,:,:) ) 660 CALL iom_put( "e3v" , e3v_n (:,:,:) ) 661 CALL iom_put( "e3w" , e3w_n (:,:,:) ) 662 CALL iom_put( "tpt_dep" , gde3w_n(:,:,:) ) 668 663 IF( iom_use("e3tdef") ) & 669 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )664 CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 670 665 671 666 ! write restart file … … 674 669 ! 675 670 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 676 671 ! 677 672 END SUBROUTINE dom_vvl_sf_swp 678 673 … … 801 796 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 802 797 ! 803 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )804 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )798 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 799 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 805 800 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 806 801 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 810 805 ! ! --------- ! 811 806 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 812 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )813 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )807 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 808 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 814 809 ! needed to restart if land processor not computed 815 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'810 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 816 811 WHERE ( tmask(:,:,:) == 0.0_wp ) 817 fse3t_n(:,:,:) = e3t_0(:,:,:)818 fse3t_b(:,:,:) = e3t_0(:,:,:)812 e3t_n(:,:,:) = e3t_0(:,:,:) 813 e3t_b(:,:,:) = e3t_0(:,:,:) 819 814 END WHERE 820 815 IF( neuler == 0 ) THEN 821 fse3t_b(:,:,:) = fse3t_n(:,:,:)816 e3t_b(:,:,:) = e3t_n(:,:,:) 822 817 ENDIF 823 818 ELSE IF( id1 > 0 ) THEN 824 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'825 IF(lwp) write(numout,*) ' fse3t_n set equal to fse3t_b.'819 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 820 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 826 821 IF(lwp) write(numout,*) 'neuler is forced to 0' 827 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )828 fse3t_n(:,:,:) = fse3t_b(:,:,:)822 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 823 e3t_n(:,:,:) = e3t_b(:,:,:) 829 824 neuler = 0 830 825 ELSE IF( id2 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'832 IF(lwp) write(numout,*) ' fse3t_b set equal to fse3t_n.'826 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 827 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 833 828 IF(lwp) write(numout,*) 'neuler is forced to 0' 834 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )835 fse3t_b(:,:,:) = fse3t_n(:,:,:)829 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 830 e3t_b(:,:,:) = e3t_n(:,:,:) 836 831 neuler = 0 837 832 ELSE 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'833 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 839 834 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 840 835 IF(lwp) write(numout,*) 'neuler is forced to 0' 841 DO jk =1,jpk842 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &843 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)&844 & + e3t_0(:,:,jk)* (1._wp -tmask(:,:,jk))836 DO jk = 1, jpk 837 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 838 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 839 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 845 840 END DO 846 fse3t_b(:,:,:) = fse3t_n(:,:,:)841 e3t_b(:,:,:) = e3t_n(:,:,:) 847 842 neuler = 0 848 843 ENDIF … … 875 870 ! 876 871 ELSE !* Initialize at "rest" 877 fse3t_b(:,:,:) = e3t_0(:,:,:)878 fse3t_n(:,:,:) = e3t_0(:,:,:)872 e3t_b(:,:,:) = e3t_0(:,:,:) 873 e3t_n(:,:,:) = e3t_0(:,:,:) 879 874 sshn(:,:) = 0.0_wp 880 875 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN … … 891 886 ! ! all cases ! 892 887 ! ! --------- ! 893 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )894 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )888 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 889 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 895 890 ! ! ----------------------- ! 896 891 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !
Note: See TracChangeset
for help on using the changeset viewer.