Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 (modified) (41 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5836 r6140 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 !!---------------------------------------------------------------------- 65 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)66 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 66 67 !! $Id$ 67 68 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 68 69 !!---------------------------------------------------------------------- 69 70 70 CONTAINS 71 71 … … 74 74 !! *** FUNCTION dom_vvl_alloc *** 75 75 !!---------------------------------------------------------------------- 76 IF( ln_vvl_zstar ) dom_vvl_alloc = 076 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 77 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 78 78 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 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 !!---------------------------------------------------------------------- 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 121 ! 122 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 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 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 181 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,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(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 193 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 194 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 195 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 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. 266 268 !!---------------------------------------------------------------------- 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 269 INTEGER, INTENT( in ) :: kt ! time step 270 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 271 ! 272 INTEGER :: ji, jj, jk ! dummy loop indices 273 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 274 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 275 LOGICAL :: ll_do_bclinic ! local logical 276 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 277 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 278 !!---------------------------------------------------------------------- 279 ! 280 IF( ln_linssh ) RETURN ! No calculation in linear free surface 281 ! 282 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 283 ! 284 CALL wrk_alloc( jpi,jpj,zht, z_scale, zwu, zwv, zhdiv ) 285 CALL wrk_alloc( jpi,jpj,jpk, ze3t ) 286 287 IF( kt == nit000 ) THEN 284 288 IF(lwp) WRITE(numout,*) 285 289 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 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------ ! 310 311 313 ! I - initialization 312 314 ! ================== … … 314 316 ! 1 - barotropic divergence 315 317 ! ------------------------- 316 zhdiv(:,:) = 0. 317 zht(:,:) = 0. 318 zhdiv(:,:) = 0._wp 319 zht(:,:) = 0._wp 318 320 DO jk = 1, jpkm1 319 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)320 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)321 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 322 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 321 323 END DO 322 324 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 325 327 ! -------------------------------------------------- 326 328 IF( ln_vvl_ztilde ) THEN 327 IF( kt .GT.nit000 ) THEN329 IF( kt > nit000 ) THEN 328 330 DO jk = 1, jpkm1 329 331 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 330 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )332 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 331 333 END DO 332 334 ENDIF 333 END 335 ENDIF 334 336 335 337 ! II - after z_tilde increments of vertical scale factors 336 338 ! ======================================================= 337 tilde_e3t_a(:,:,:) = 0. 0_wp ! tilde_e3t_a used to store tendency terms339 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 338 340 339 341 ! 1 - High frequency divergence term … … 341 343 IF( ln_vvl_ztilde ) THEN ! z_tilde case 342 344 DO jk = 1, jpkm1 343 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )345 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 344 346 END DO 345 347 ELSE ! layer case 346 348 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 349 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 350 END DO 351 ENDIF 350 352 351 353 ! 2 - Restoring term (z-tilde case only) … … 355 357 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 356 358 END DO 357 END 359 ENDIF 358 360 359 361 ! 3 - Thickness diffusion term 360 362 ! ---------------------------- 361 zwu(:,:) = 0.0_wp 362 zwv(:,:) = 0.0_wp 363 ! a - first derivative: diffusive fluxes 364 DO jk = 1, jpkm1 363 zwu(:,:) = 0._wp 364 zwv(:,:) = 0._wp 365 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 365 366 DO jj = 1, jpjm1 366 367 DO ji = 1, fs_jpim1 ! vector opt. … … 374 375 END DO 375 376 END DO 376 ! b - correction for last oceanic u-v points 377 DO jj = 1, jpj 377 DO jj = 1, jpj ! b - correction for last oceanic u-v points 378 378 DO ji = 1, jpi 379 379 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 381 381 END DO 382 382 END DO 383 ! c - second derivative: divergence of diffusive fluxes 384 DO jk = 1, jpkm1 383 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 385 384 DO jj = 2, jpjm1 386 385 DO ji = fs_2, fs_jpim1 ! vector opt. … … 391 390 END DO 392 391 END DO 393 ! d - thickness diffusion transport: boundary conditions394 ! (stored for tracer advction and continuity equation)392 ! ! d - thickness diffusion transport: boundary conditions 393 ! (stored for tracer advction and continuity equation) 395 394 CALL lbc_lnk( un_td , 'U' , -1._wp) 396 395 CALL lbc_lnk( vn_td , 'V' , -1._wp) … … 410 409 ! Maximum deformation control 411 410 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 412 ze3t(:,:,jpk) = 0. 0_wp411 ze3t(:,:,jpk) = 0._wp 413 412 DO jk = 1, jpkm1 414 413 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 419 418 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 420 419 ! - 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 ) ) THEN420 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 422 421 IF( lk_mpp ) THEN 423 422 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 462 461 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 463 462 DO jk = 1, jpkm1 464 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)463 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 465 464 END DO 466 465 … … 470 469 ! ! ---baroclinic part--------- ! 471 470 DO jk = 1, jpkm1 472 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)471 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 473 472 END DO 474 473 ENDIF … … 485 484 zht(:,:) = 0.0_wp 486 485 DO jk = 1, jpkm1 487 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)486 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 488 487 END DO 489 488 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 490 489 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_tmax490 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 492 491 ! 493 492 zht(:,:) = 0.0_wp 494 493 DO jk = 1, jpkm1 495 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)494 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 496 495 END DO 497 496 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 498 497 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_tmax498 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 500 499 ! 501 500 zht(:,:) = 0.0_wp 502 501 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)502 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 504 503 END DO 505 504 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 506 505 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_tmax506 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 508 507 ! 509 508 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 524 523 ! *********************************** ! 525 524 526 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )527 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )525 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 526 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 528 527 529 528 ! *********************************** ! … … 531 530 ! *********************************** ! 532 531 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)532 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 533 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 534 DO jk = 2, jpkm1 535 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 536 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 538 537 END DO 539 538 ! ! 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 539 !!gm BUG ? don't understand the use of umask_i here ..... 540 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 541 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 542 ! 543 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) 544 CALL wrk_dealloc( jpi,jpj,jpk, ze3t ) 545 ! 546 546 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 547 547 ! 548 548 END SUBROUTINE dom_vvl_sf_nxt 549 549 … … 561 561 !! - recompute depths and water height fields 562 562 !! 563 !! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step563 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 564 564 !! - Recompute: 565 !! fse3(u/v)_b566 !! fse3w_n567 !! fse3(u/v)w_b568 !! fse3(u/v)w_n569 !! fsdept_n, fsdepw_n and fsde3w_n565 !! e3(u/v)_b 566 !! e3w_n 567 !! e3(u/v)w_b 568 !! e3(u/v)w_n 569 !! gdept_n, gdepw_n and gde3w_n 570 570 !! h(u/v) and h(u/v)r 571 571 !! … … 573 573 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 574 574 !!---------------------------------------------------------------------- 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 575 INTEGER, INTENT( in ) :: kt ! time step 576 ! 577 INTEGER :: ji, jj, jk ! dummy loop indices 578 REAL(wp) :: zcoef ! local scalar 579 !!---------------------------------------------------------------------- 580 ! 581 IF( ln_linssh ) RETURN ! No calculation in linear free surface 582 ! 582 583 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 583 584 ! … … 590 591 ! Time filter and swap of scale factors 591 592 ! ===================================== 592 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.593 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 593 594 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 594 595 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 600 601 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 601 602 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(:,:,:)603 gdept_b(:,:,:) = gdept_n(:,:,:) 604 gdepw_b(:,:,:) = gdepw_n(:,:,:) 605 606 e3t_n(:,:,:) = e3t_a(:,:,:) 607 e3u_n(:,:,:) = e3u_a(:,:,:) 608 e3v_n(:,:,:) = e3v_a(:,:,:) 608 609 609 610 ! Compute all missing vertical scale factor and depths … … 611 612 ! Horizontal scale factor interpolations 612 613 ! -------------------------------------- 613 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt614 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 614 615 ! - JC - hu_b, hv_b, hur_b, hvr_b also 615 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 616 617 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 618 616 619 ! 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 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' ) 626 627 ! t- and w- points depth (set the isf depth as it is in the initial step) 628 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 629 gdepw_n(:,:,1) = 0.0_wp 630 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 631 631 DO jk = 2, jpk 632 632 DO jj = 1,jpj … … 635 635 ! 1 for jk = mikt 636 636 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)637 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 638 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 639 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 640 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 641 641 END DO 642 642 END DO 643 643 END DO 644 644 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) 645 ! Local depth and Inverse of the local depth of the water 646 ! ------------------------------------------------------- 647 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 648 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 649 ! 650 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 651 DO jk = 2, jpkm1 652 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 659 653 END DO 660 654 661 655 ! Write outputs 662 656 ! ============= 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(:,:,:) )657 CALL iom_put( "e3t", e3t_n(:,:,:) ) 658 CALL iom_put( "e3u", e3u_n(:,:,:) ) 659 CALL iom_put( "e3v", e3v_n(:,:,:) ) 660 CALL iom_put( "e3w", e3w_n(:,:,:) ) 661 CALL iom_put( "tpt_dep", gde3w_n(:,:,:) ) 668 662 IF( iom_use("e3tdef") ) & 669 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100* tmask(:,:,:) ) ** 2 )663 CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 ) 670 664 671 665 ! write restart file 672 666 ! ================== 673 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )674 ! 675 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')676 667 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 668 ! 669 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 670 ! 677 671 END SUBROUTINE dom_vvl_sf_swp 678 672 … … 696 690 !!---------------------------------------------------------------------- 697 691 ! 698 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol')692 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 699 693 ! 700 694 SELECT CASE ( pout ) !== type of interpolation ==! … … 770 764 END SELECT 771 765 ! 772 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol')766 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 773 767 ! 774 768 END SUBROUTINE dom_vvl_interpol … … 801 795 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 802 796 ! 803 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )804 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )797 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 798 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 805 799 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 806 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 810 804 ! ! --------- ! 811 805 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(:,:,:) )806 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 807 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 814 808 ! 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'809 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 816 810 WHERE ( tmask(:,:,:) == 0.0_wp ) 817 fse3t_n(:,:,:) = e3t_0(:,:,:)818 fse3t_b(:,:,:) = e3t_0(:,:,:)811 e3t_n(:,:,:) = e3t_0(:,:,:) 812 e3t_b(:,:,:) = e3t_0(:,:,:) 819 813 END WHERE 820 814 IF( neuler == 0 ) THEN 821 fse3t_b(:,:,:) = fse3t_n(:,:,:)815 e3t_b(:,:,:) = e3t_n(:,:,:) 822 816 ENDIF 823 817 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.'818 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 819 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 826 820 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(:,:,:)821 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 822 e3t_n(:,:,:) = e3t_b(:,:,:) 829 823 neuler = 0 830 824 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.'825 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 826 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 833 827 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(:,:,:)828 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 829 e3t_b(:,:,:) = e3t_n(:,:,:) 836 830 neuler = 0 837 831 ELSE 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'832 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 839 833 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 840 834 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))835 DO jk = 1, jpk 836 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 837 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 838 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 845 839 END DO 846 fse3t_b(:,:,:) = fse3t_n(:,:,:)840 e3t_b(:,:,:) = e3t_n(:,:,:) 847 841 neuler = 0 848 842 ENDIF … … 875 869 ! 876 870 ELSE !* Initialize at "rest" 877 fse3t_b(:,:,:) = e3t_0(:,:,:)878 fse3t_n(:,:,:) = e3t_0(:,:,:)871 e3t_b(:,:,:) = e3t_0(:,:,:) 872 e3t_n(:,:,:) = e3t_0(:,:,:) 879 873 sshn(:,:) = 0.0_wp 880 874 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN … … 891 885 ! ! all cases ! 892 886 ! ! --------- ! 893 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )894 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )887 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 888 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 895 889 ! ! ----------------------- ! 896 890 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 975 969 ! 976 970 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' )971 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 978 972 ! 979 973 IF(lwp) THEN ! Print the choice … … 989 983 ! 990 984 #if defined key_agrif 991 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )985 IF(.NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 992 986 #endif 993 987 ! … … 996 990 !!====================================================================== 997 991 END MODULE domvvl 998 999 1000
Note: See TracChangeset
for help on using the changeset viewer.
