- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5836 r7351 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 wet_dry ! wetting and drying 26 USE restart ! ocean restart 27 ! 24 28 USE in_out_manager ! I/O manager 25 29 USE iom ! I/O manager library 26 USE restart ! ocean restart27 30 USE lib_mpp ! distributed memory computing library 28 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 60 63 61 64 !! * Substitutions 62 # include "domzgr_substitute.h90"63 65 # include "vectopt_loop_substitute.h90" 64 66 !!---------------------------------------------------------------------- 65 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)67 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 66 68 !! $Id$ 67 69 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 68 70 !!---------------------------------------------------------------------- 69 70 71 CONTAINS 71 72 … … 74 75 !! *** FUNCTION dom_vvl_alloc *** 75 76 !!---------------------------------------------------------------------- 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 077 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 78 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 79 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 81 82 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 82 83 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_wp84 un_td = 0._wp 85 vn_td = 0._wp 85 86 ENDIF 86 87 IF( ln_vvl_ztilde ) THEN … … 89 90 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 90 91 ENDIF 91 92 ! 92 93 END FUNCTION dom_vvl_alloc 93 94 … … 103 104 !! - interpolate scale factors 104 105 !! 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_n106 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 107 !! - Regrid: e3(u/v)_n 108 !! e3(u/v)_b 109 !! e3w_n 110 !! e3(u/v)w_b 111 !! e3(u/v)w_n 112 !! gdept_n, gdepw_n and gde3w_n 112 113 !! - h(t/u/v)_0 113 114 !! - frq_rst_e3t and frq_rst_hdv … … 115 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 117 !!---------------------------------------------------------------------- 117 USE phycst, ONLY : rpi, rsmall, rad 118 !! * Local declarations 119 INTEGER :: ji,jj,jk 118 INTEGER :: ji, jj, jk 120 119 INTEGER :: ii0, ii1, ij0, ij1 121 120 REAL(wp):: zcoef 122 121 !!---------------------------------------------------------------------- 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 122 ! 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 ! 125 125 IF(lwp) WRITE(numout,*) 126 126 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 127 127 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 ! ====================== 128 ! 129 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! 131 ! ! Allocate module arrays 135 132 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 ! ============================================================================ 133 ! 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 139 135 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 136 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 137 ! 138 ! !== Set of all other vertical scale factors ==! (now and before) 139 ! ! Horizontal interpolation of e3t 140 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 141 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 142 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 143 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 144 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 145 ! ! Vertical interpolation of e3t,u,v 146 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 147 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 148 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 149 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 150 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 151 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 152 ! 153 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 154 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 155 gdepw_n(:,:,1) = 0.0_wp 156 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 157 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 158 gdepw_b(:,:,1) = 0.0_wp 159 DO jk = 2, jpk ! vertical sum 169 160 DO jj = 1,jpj 170 161 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)) 162 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 163 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 164 ! ! 0.5 where jk = mikt 165 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 166 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 167 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 168 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 169 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 170 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 171 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 172 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 173 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 182 174 END DO 183 175 END DO 184 176 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) 177 ! 178 ! !== thickness of the water column !! (ocean portion only) 179 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 180 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 181 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 182 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 183 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 184 DO jk = 2, jpkm1 185 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 186 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 187 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 188 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 189 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 193 190 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 ! ============================================ 191 ! 192 ! !== inverse of water column thickness ==! (u- and v- points) 193 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 194 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 195 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 196 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 197 198 ! !== z_tilde coordinate case ==! (Restoring frequencies) 199 199 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 200 !!gm : idea: add here a READ in a file of custumized restoring frequency 201 ! ! Values in days provided via the namelist 202 ! ! use rsmall to avoid possible division by zero errors with faulty settings 203 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 204 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 205 ! 206 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 207 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 208 frq_rst_hdv(:,:) = 1._wp / rdt 207 209 ENDIF 208 IF ( ln_vvl_zstar_at_eqtor ) THEN 210 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 209 211 DO jj = 1, jpj 210 212 DO ji = 1, jpi 213 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 211 214 IF( ABS(gphit(ji,jj)) >= 6.) THEN 212 215 ! values outside the equatorial band and transition zone (ztilde) 213 216 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 214 217 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 218 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 216 219 ! values inside the equatorial band (ztilde as zstar) 217 220 frq_rst_e3t(ji,jj) = 0.0_wp 218 221 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)222 ELSE ! transition band (2.5 to 6 degrees N/S) 223 ! ! (linearly transition from z-tilde to z-star) 221 224 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 222 225 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 229 232 END DO 230 233 END DO 231 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 232 ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2234 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 235 ii0 = 103 ; ii1 = 111 233 236 ij0 = 128 ; ij1 = 135 ; 234 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 237 240 ENDIF 238 241 ENDIF 239 242 ! 240 243 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 241 244 ! 242 245 END SUBROUTINE dom_vvl_init 243 246 … … 261 264 !! - tilde_e3t_a: after increment of vertical scale factor 262 265 !! in z_tilde case 263 !! - fse3(t/u/v)_a266 !! - e3(t/u/v)_a 264 267 !! 265 268 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 266 269 !!---------------------------------------------------------------------- 267 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 268 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 269 !! * Arguments 270 INTEGER, INTENT( in ) :: kt ! time step 271 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 272 !! * Local declarations 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 275 REAL(wp) :: z2dt ! temporary scalars 276 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 277 LOGICAL :: ll_do_bclinic ! temporary logical 278 !!---------------------------------------------------------------------- 279 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 ) 282 283 IF(kt == nit000) THEN 270 INTEGER, INTENT( in ) :: kt ! time step 271 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 272 ! 273 INTEGER :: ji, jj, jk ! dummy loop indices 274 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 275 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 276 LOGICAL :: ll_do_bclinic ! local logical 277 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 278 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 279 !!---------------------------------------------------------------------- 280 ! 281 IF( ln_linssh ) RETURN ! No calculation in linear free surface 282 ! 283 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 284 ! 285 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv ) 286 CALL wrk_alloc( jpi,jpj,jpk, ze3t ) 287 288 IF( kt == nit000 ) THEN 284 289 IF(lwp) WRITE(numout,*) 285 290 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 289 294 ll_do_bclinic = .TRUE. 290 295 IF( PRESENT(kcall) ) THEN 291 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.296 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 292 297 ENDIF 293 298 … … 295 300 ! After acale factors at t-points ! 296 301 ! ******************************* ! 297 298 302 ! ! --------------------------------------------- ! 299 303 ! ! z_star coordinate and barotropic z-tilde part ! 300 304 ! ! --------------------------------------------- ! 301 305 ! 302 306 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 303 307 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)308 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 309 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 306 310 END DO 307 311 ! 308 312 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 309 313 ! ! ------baroclinic part------ ! 310 311 314 ! I - initialization 312 315 ! ================== … … 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(:,:) … … 419 419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 420 420 ! - ML - test: for the moment, stop simulation for too large e3_t variations 421 IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT.- rn_zdef_max ) ) THEN421 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 422 422 IF( lk_mpp ) THEN 423 423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 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(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 542 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 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 580 !!---------------------------------------------------------------------- 581 576 INTEGER, INTENT( in ) :: kt ! time step 577 ! 578 INTEGER :: ji, jj, jk ! dummy loop indices 579 REAL(wp) :: zcoef ! local scalar 580 !!---------------------------------------------------------------------- 581 ! 582 IF( ln_linssh ) RETURN ! No calculation in linear free surface 583 ! 582 584 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 583 585 ! … … 590 592 ! Time filter and swap of scale factors 591 593 ! ===================================== 592 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.594 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 593 595 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 594 596 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 600 602 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 601 603 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(:,:,:)604 gdept_b(:,:,:) = gdept_n(:,:,:) 605 gdepw_b(:,:,:) = gdepw_n(:,:,:) 606 607 e3t_n(:,:,:) = e3t_a(:,:,:) 608 e3u_n(:,:,:) = e3u_a(:,:,:) 609 e3v_n(:,:,:) = e3v_a(:,:,:) 608 610 609 611 ! Compute all missing vertical scale factor and depths … … 611 613 ! Horizontal scale factor interpolations 612 614 ! -------------------------------------- 613 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt615 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 614 616 ! - JC - hu_b, hv_b, hur_b, hvr_b also 615 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 617 618 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 619 616 620 ! Vertical scale factor interpolations 617 ! ------------------------------------ 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' ) 624 ! t- and w- points depth 625 ! ---------------------- 626 ! 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_wp 629 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 630 621 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 622 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 623 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 624 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 625 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 626 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 627 628 ! t- and w- points depth (set the isf depth as it is in the initial step) 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(:,:) 631 632 DO jk = 2, jpk 632 633 DO jj = 1,jpj … … 635 636 ! 1 for jk = mikt 636 637 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)638 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 639 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 640 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 641 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 641 642 END DO 642 643 END DO 643 644 END DO 644 645 645 ! Local depth and Inverse of the local depth of the water column at u- and v- points 646 ! ---------------------------------------------------------------------------------- 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) 646 ! Local depth and Inverse of the local depth of the water 647 ! ------------------------------------------------------- 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 661 ! Write outputs662 ! =============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 (:,:,:) )668 IF( iom_use("e3tdef") ) &669 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )670 655 671 656 ! write restart file 672 657 ! ================== 673 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )674 ! 675 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')676 658 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 659 ! 660 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 661 ! 677 662 END SUBROUTINE dom_vvl_sf_swp 678 663 … … 693 678 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 694 679 ! 695 INTEGER :: ji, jj, jk ! dummy loop indices 696 !!---------------------------------------------------------------------- 697 ! 698 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 680 INTEGER :: ji, jj, jk ! dummy loop indices 681 REAL(wp) :: zlnwd ! =1./0. when ln_wd = T/F 682 !!---------------------------------------------------------------------- 683 ! 684 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 685 ! 686 IF(ln_wd) THEN 687 zlnwd = 1.0_wp 688 ELSE 689 zlnwd = 0.0_wp 690 END IF 699 691 ! 700 692 SELECT CASE ( pout ) !== type of interpolation ==! … … 704 696 DO jj = 1, jpjm1 705 697 DO ji = 1, fs_jpim1 ! vector opt. 706 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj)&698 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 707 699 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 708 700 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) … … 717 709 DO jj = 1, jpjm1 718 710 DO ji = 1, fs_jpim1 ! vector opt. 719 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj)&711 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 720 712 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 721 713 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) … … 730 722 DO jj = 1, jpjm1 731 723 DO ji = 1, fs_jpim1 ! vector opt. 732 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1e2f(ji,jj) & 724 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 725 & * r1_e1e2f(ji,jj) & 733 726 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 734 727 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) … … 745 738 !!gm BUG? use here wmask in case of ISF ? to be checked 746 739 DO jk = 2, jpk 747 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 748 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 740 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 741 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 742 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 743 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 749 744 END DO 750 745 ! … … 755 750 !!gm BUG? use here wumask in case of ISF ? to be checked 756 751 DO jk = 2, jpk 757 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 758 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 752 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 753 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 754 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 755 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 759 756 END DO 760 757 ! … … 765 762 !!gm BUG? use here wvmask in case of ISF ? to be checked 766 763 DO jk = 2, jpk 767 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 768 & + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 764 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 765 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 766 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 767 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 769 768 END DO 770 769 END SELECT 771 770 ! 772 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')771 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 773 772 ! 774 773 END SUBROUTINE dom_vvl_interpol … … 790 789 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 791 790 ! 792 INTEGER :: j k791 INTEGER :: ji, jj, jk 793 792 INTEGER :: id1, id2, id3, id4, id5 ! local integers 794 793 !!---------------------------------------------------------------------- … … 801 800 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 802 801 ! 803 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )804 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )802 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 803 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 805 804 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 806 805 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 810 809 ! ! --------- ! 811 810 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(:,:,:) )811 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 812 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 814 813 ! 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'814 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 816 815 WHERE ( tmask(:,:,:) == 0.0_wp ) 817 fse3t_n(:,:,:) = e3t_0(:,:,:)818 fse3t_b(:,:,:) = e3t_0(:,:,:)816 e3t_n(:,:,:) = e3t_0(:,:,:) 817 e3t_b(:,:,:) = e3t_0(:,:,:) 819 818 END WHERE 820 819 IF( neuler == 0 ) THEN 821 fse3t_b(:,:,:) = fse3t_n(:,:,:)820 e3t_b(:,:,:) = e3t_n(:,:,:) 822 821 ENDIF 823 822 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.'823 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 824 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 826 825 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(:,:,:)826 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 827 e3t_n(:,:,:) = e3t_b(:,:,:) 829 828 neuler = 0 830 829 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.'830 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 831 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 833 832 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(:,:,:)833 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 834 e3t_b(:,:,:) = e3t_n(:,:,:) 836 835 neuler = 0 837 836 ELSE 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'837 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 839 838 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 840 839 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))840 DO jk = 1, jpk 841 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 842 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 843 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 845 844 END DO 846 fse3t_b(:,:,:) = fse3t_n(:,:,:)845 e3t_b(:,:,:) = e3t_n(:,:,:) 847 846 neuler = 0 848 847 ENDIF … … 875 874 ! 876 875 ELSE !* Initialize at "rest" 877 fse3t_b(:,:,:) = e3t_0(:,:,:)878 fse3t_n(:,:,:) = e3t_0(:,:,:)876 e3t_b(:,:,:) = e3t_0(:,:,:) 877 e3t_n(:,:,:) = e3t_0(:,:,:) 879 878 sshn(:,:) = 0.0_wp 879 880 IF( ln_wd ) THEN 881 DO jj = 1, jpj 882 DO ji = 1, jpi 883 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 884 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 885 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 886 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 887 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 888 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 889 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 890 ENDIF 891 ENDDO 892 ENDDO 893 END IF 894 880 895 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 881 896 tilde_e3t_b(:,:,:) = 0.0_wp … … 891 906 ! ! all cases ! 892 907 ! ! --------- ! 893 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )894 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )908 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 909 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 895 910 ! ! ----------------------- ! 896 911 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 975 990 ! 976 991 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 977 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )992 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 978 993 ! 979 994 IF(lwp) THEN ! Print the choice … … 989 1004 ! 990 1005 #if defined key_agrif 991 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )1006 IF(.NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 992 1007 #endif 993 1008 ! … … 996 1011 !!====================================================================== 997 1012 END MODULE domvvl 998 999 1000
Note: See TracChangeset
for help on using the changeset viewer.
