- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4624 r6225 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 9 !! vvl option includes z_star and z_tilde coordinates 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 11 !!---------------------------------------------------------------------- 11 !! 'key_vvl' variable volume 12 !!---------------------------------------------------------------------- 12 13 13 !!---------------------------------------------------------------------- 14 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness … … 18 18 !! dom_vvl_rst : read/write restart file 19 19 !! dom_vvl_ctl : Check the vvl options 20 !! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors21 !! : to account for manual changes to e[1,2][u,v] in some Straits22 20 !!---------------------------------------------------------------------- 23 !! * Modules used24 21 USE oce ! ocean dynamics and tracers 22 USE phycst ! physical constant 25 23 USE dom_oce ! ocean space and time domain 26 24 USE sbc_oce ! ocean surface boundary condition 25 USE wet_dry ! wetting and drying 26 USE restart ! ocean restart 27 ! 27 28 USE in_out_manager ! I/O manager 28 29 USE iom ! I/O manager library 29 USE restart ! ocean restart30 30 USE lib_mpp ! distributed memory computing library 31 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 36 36 PRIVATE 37 37 38 !! * Routine accessibility39 38 PUBLIC dom_vvl_init ! called by domain.F90 40 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 41 40 PUBLIC dom_vvl_sf_swp ! called by step.F90 42 41 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 43 PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol 44 45 !!* Namelist nam_vvl 46 LOGICAL , PUBLIC :: ln_vvl_zstar ! zstar vertical coordinate 47 LOGICAL , PUBLIC :: ln_vvl_ztilde ! ztilde vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_layer ! level vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar ! ztilde vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_kepe ! kinetic/potential energy transfer 52 ! ! conservation: not used yet 53 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 54 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 55 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 56 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 57 LOGICAL , PUBLIC :: ln_vvl_dbg ! debug control prints 58 59 !! * Module variables 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 64 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 42 43 ! !!* Namelist nam_vvl 44 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 45 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 46 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 47 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 50 ! ! conservation: not used yet 51 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 52 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 53 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 54 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 55 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 56 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 58 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 62 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 66 63 67 64 !! * Substitutions 68 # include "domzgr_substitute.h90"69 65 # include "vectopt_loop_substitute.h90" 70 66 !!---------------------------------------------------------------------- 71 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)67 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 72 68 !! $Id$ 73 69 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 74 70 !!---------------------------------------------------------------------- 75 76 71 CONTAINS 77 72 … … 80 75 !! *** FUNCTION dom_vvl_alloc *** 81 76 !!---------------------------------------------------------------------- 82 IF( ln_vvl_zstar ) dom_vvl_alloc = 077 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 83 78 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 84 79 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 87 82 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 88 83 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 89 un_td = 0. 0_wp90 vn_td = 0. 0_wp84 un_td = 0._wp 85 vn_td = 0._wp 91 86 ENDIF 92 87 IF( ln_vvl_ztilde ) THEN … … 95 90 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 96 91 ENDIF 97 92 ! 98 93 END FUNCTION dom_vvl_alloc 99 94 … … 109 104 !! - interpolate scale factors 110 105 !! 111 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)112 !! - Regrid: fse3(u/v)_n113 !! fse3(u/v)_b114 !! fse3w_n115 !! fse3(u/v)w_b116 !! fse3(u/v)w_n117 !! 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 118 113 !! - h(t/u/v)_0 119 114 !! - frq_rst_e3t and frq_rst_hdv … … 121 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 122 117 !!---------------------------------------------------------------------- 123 USE phycst, ONLY : rpi, rsmall, rad 124 !! * Local declarations 125 INTEGER :: ji,jj,jk 118 INTEGER :: ji, jj, jk 126 119 INTEGER :: ii0, ii1, ij0, ij1 127 !!---------------------------------------------------------------------- 128 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 129 120 REAL(wp):: zcoef 121 !!---------------------------------------------------------------------- 122 ! 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 ! 130 125 IF(lwp) WRITE(numout,*) 131 126 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 132 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 133 134 ! choose vertical coordinate (z_star, z_tilde or layer) 135 ! ========================== 136 CALL dom_vvl_ctl 137 138 ! Allocate module arrays 139 ! ====================== 128 ! 129 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! 131 ! ! Allocate module arrays 140 132 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 141 142 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 143 ! ============================================================================ 133 ! 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 144 135 CALL dom_vvl_rst( nit000, 'READ' ) 145 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 146 147 ! Reconstruction of all vertical scale factors at now and before time steps 148 ! ============================================================================= 149 ! Horizontal scale factor interpolations 150 ! -------------------------------------- 151 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 152 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 153 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 154 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 155 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 156 ! Vertical scale factor interpolations 157 ! ------------------------------------ 158 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 159 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 160 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 161 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 162 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 163 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 164 ! t- and w- points depth 165 ! ---------------------- 166 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 167 fsdepw_n(:,:,1) = 0.0_wp 168 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 169 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 170 fsdepw_b(:,:,1) = 0.0_wp 171 DO jk = 2, jpk 172 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 173 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 174 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 175 fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk) 176 fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1) 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 160 DO jj = 1,jpj 161 DO ji = 1,jpi 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)) 174 END DO 175 END DO 177 176 END DO 178 179 ! Before depth and Inverse of the local depth of the water column at u- and v- points 180 ! ----------------------------------------------------------------------------------- 181 hu_b(:,:) = 0. 182 hv_b(:,:) = 0. 183 DO jk = 1, jpkm1 184 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 185 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) 186 190 END DO 187 hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) ) 188 hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) ) 189 190 ! Restoring frequencies for z_tilde coordinate 191 ! ============================================ 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) 192 199 IF( ln_vvl_ztilde ) THEN 193 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 194 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 195 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 196 IF( ln_vvl_ztilde_as_zstar ) THEN 197 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 198 frq_rst_e3t(:,:) = 0.0_wp 199 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 200 209 ENDIF 201 IF ( ln_vvl_zstar_at_eqtor ) THEN 210 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 202 211 DO jj = 1, jpj 203 212 DO ji = 1, jpi 213 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 204 214 IF( ABS(gphit(ji,jj)) >= 6.) THEN 205 215 ! values outside the equatorial band and transition zone (ztilde) 206 216 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 207 217 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 208 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 218 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 209 219 ! values inside the equatorial band (ztilde as zstar) 210 220 frq_rst_e3t(ji,jj) = 0.0_wp 211 221 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 212 ELSE 213 ! 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) 214 224 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 215 225 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 222 232 END DO 223 233 END DO 224 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 225 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 226 236 ij0 = 128 ; ij1 = 135 ; 227 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 230 240 ENDIF 231 241 ENDIF 232 242 ! 233 243 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 234 244 ! 235 245 END SUBROUTINE dom_vvl_init 236 246 … … 254 264 !! - tilde_e3t_a: after increment of vertical scale factor 255 265 !! in z_tilde case 256 !! - fse3(t/u/v)_a266 !! - e3(t/u/v)_a 257 267 !! 258 268 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 259 269 !!---------------------------------------------------------------------- 260 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 261 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 262 !! * Arguments 263 INTEGER, INTENT( in ) :: kt ! time step 264 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 265 !! * Local declarations 266 INTEGER :: ji, jj, jk ! dummy loop indices 267 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 268 REAL(wp) :: z2dt ! temporary scalars 269 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 270 LOGICAL :: ll_do_bclinic ! temporary logical 271 !!---------------------------------------------------------------------- 272 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 273 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 274 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 275 276 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 277 289 IF(lwp) WRITE(numout,*) 278 290 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 282 294 ll_do_bclinic = .TRUE. 283 295 IF( PRESENT(kcall) ) THEN 284 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.296 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 285 297 ENDIF 286 298 … … 288 300 ! After acale factors at t-points ! 289 301 ! ******************************* ! 290 291 302 ! ! --------------------------------------------- ! 292 303 ! ! z_star coordinate and barotropic z-tilde part ! 293 304 ! ! --------------------------------------------- ! 294 295 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )305 ! 306 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 296 307 DO jk = 1, jpkm1 297 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)298 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) 299 310 END DO 300 311 ! 301 312 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 302 313 ! ! ------baroclinic part------ ! 303 304 314 ! I - initialization 305 315 ! ================== … … 307 317 ! 1 - barotropic divergence 308 318 ! ------------------------- 309 zhdiv(:,:) = 0. 310 zht(:,:) = 0. 319 zhdiv(:,:) = 0._wp 320 zht(:,:) = 0._wp 311 321 DO jk = 1, jpkm1 312 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)313 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)314 END DO 315 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask (:,:,1) )322 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 323 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 324 END DO 325 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 316 326 317 327 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 318 328 ! -------------------------------------------------- 319 329 IF( ln_vvl_ztilde ) THEN 320 IF( kt .GT.nit000 ) THEN330 IF( kt > nit000 ) THEN 321 331 DO jk = 1, jpkm1 322 332 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 323 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )333 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 324 334 END DO 325 335 ENDIF 326 END 336 ENDIF 327 337 328 338 ! II - after z_tilde increments of vertical scale factors 329 339 ! ======================================================= 330 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 331 341 332 342 ! 1 - High frequency divergence term … … 334 344 IF( ln_vvl_ztilde ) THEN ! z_tilde case 335 345 DO jk = 1, jpkm1 336 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) ) 337 347 END DO 338 348 ELSE ! layer case 339 349 DO jk = 1, jpkm1 340 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:))341 END DO 342 END 350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 351 END DO 352 ENDIF 343 353 344 354 ! 2 - Restoring term (z-tilde case only) … … 348 358 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 349 359 END DO 350 END 360 ENDIF 351 361 352 362 ! 3 - Thickness diffusion term 353 363 ! ---------------------------- 354 zwu(:,:) = 0.0_wp 355 zwv(:,:) = 0.0_wp 356 ! a - first derivative: diffusive fluxes 357 DO jk = 1, jpkm1 364 zwu(:,:) = 0._wp 365 zwv(:,:) = 0._wp 366 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 358 367 DO jj = 1, jpjm1 359 368 DO ji = 1, fs_jpim1 ! vector opt. 360 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj)&361 &* ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )362 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj)&363 &* ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) )369 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 370 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 371 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 372 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 364 373 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 365 374 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 367 376 END DO 368 377 END DO 369 ! b - correction for last oceanic u-v points 370 DO jj = 1, jpj 378 DO jj = 1, jpj ! b - correction for last oceanic u-v points 371 379 DO ji = 1, jpi 372 380 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 374 382 END DO 375 383 END DO 376 ! c - second derivative: divergence of diffusive fluxes 377 DO jk = 1, jpkm1 384 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 378 385 DO jj = 2, jpjm1 379 386 DO ji = fs_2, fs_jpim1 ! vector opt. 380 387 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 381 388 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 382 & ) * r1_e1 2t(ji,jj)389 & ) * r1_e1e2t(ji,jj) 383 390 END DO 384 391 END DO 385 392 END DO 386 ! d - thickness diffusion transport: boundary conditions387 ! (stored for tracer advction and continuity equation)388 CALL lbc_lnk( un_td , 'U' , -1. )389 CALL lbc_lnk( vn_td , 'V' , -1. )393 ! ! d - thickness diffusion transport: boundary conditions 394 ! (stored for tracer advction and continuity equation) 395 CALL lbc_lnk( un_td , 'U' , -1._wp) 396 CALL lbc_lnk( vn_td , 'V' , -1._wp) 390 397 391 398 ! 4 - Time stepping of baroclinic scale factors … … 398 405 z2dt = 2.0_wp * rdt 399 406 ENDIF 400 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )407 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 401 408 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 402 409 403 410 ! Maximum deformation control 404 411 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 405 ze3t(:,:,jpk) = 0. 0_wp412 ze3t(:,:,jpk) = 0._wp 406 413 DO jk = 1, jpkm1 407 414 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 412 419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 413 420 ! - ML - test: for the moment, stop simulation for too large e3_t variations 414 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 415 422 IF( lk_mpp ) THEN 416 423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 453 460 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 454 461 END DO 455 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )462 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 456 463 DO jk = 1, jpkm1 457 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) 458 465 END DO 459 466 … … 463 470 ! ! ---baroclinic part--------- ! 464 471 DO jk = 1, jpkm1 465 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)472 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 466 473 END DO 467 474 ENDIF … … 478 485 zht(:,:) = 0.0_wp 479 486 DO jk = 1, jpkm1 480 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)487 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 481 488 END DO 482 489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 483 490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 484 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 485 492 ! 486 493 zht(:,:) = 0.0_wp 487 494 DO jk = 1, jpkm1 488 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)495 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 489 496 END DO 490 497 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 491 498 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 492 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 493 500 ! 494 501 zht(:,:) = 0.0_wp 495 502 DO jk = 1, jpkm1 496 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)503 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 497 504 END DO 498 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 499 506 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 500 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 501 508 ! 502 509 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 517 524 ! *********************************** ! 518 525 519 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )520 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' ) 521 528 522 529 ! *********************************** ! … … 524 531 ! *********************************** ! 525 532 526 hu_a(:,:) = 0._wp ! Ocean depth at U-points527 hv_a(:,:) = 0._wp ! Ocean depth at V-points528 DO jk = 1, jpkm1529 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)530 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) 531 538 END DO 532 539 ! ! Inverse of the local depth 533 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 534 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 535 536 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 537 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 538 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 ! 539 547 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 540 548 ! 541 549 END SUBROUTINE dom_vvl_sf_nxt 542 550 … … 554 562 !! - recompute depths and water height fields 555 563 !! 556 !! ** 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 557 565 !! - Recompute: 558 !! fse3(u/v)_b559 !! fse3w_n560 !! fse3(u/v)w_b561 !! fse3(u/v)w_n562 !! 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 563 571 !! h(u/v) and h(u/v)r 564 572 !! … … 566 574 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 567 575 !!---------------------------------------------------------------------- 568 !! * Arguments 569 INTEGER, INTENT( in ) :: kt ! time step 570 !! * Local declarations 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 572 INTEGER :: jk ! dummy loop indices 573 !!---------------------------------------------------------------------- 574 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 ! 575 584 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 576 !577 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def )578 585 ! 579 586 IF( kt == nit000 ) THEN … … 585 592 ! Time filter and swap of scale factors 586 593 ! ===================================== 587 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.594 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 588 595 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 589 596 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 595 602 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 596 603 ENDIF 597 fsdept_b(:,:,:) = fsdept_n(:,:,:)598 fsdepw_b(:,:,:) = fsdepw_n(:,:,:)599 600 fse3t_n(:,:,:) = fse3t_a(:,:,:)601 fse3u_n(:,:,:) = fse3u_a(:,:,:)602 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(:,:,:) 603 610 604 611 ! Compute all missing vertical scale factor and depths … … 606 613 ! Horizontal scale factor interpolations 607 614 ! -------------------------------------- 608 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt615 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 609 616 ! - JC - hu_b, hv_b, hur_b, hvr_b also 610 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 617 618 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 619 611 620 ! Vertical scale factor interpolations 612 ! ------------------------------------ 613 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 614 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 615 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 616 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 617 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 618 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 619 ! t- and w- points depth 620 ! ---------------------- 621 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 622 fsdepw_n(:,:,1) = 0.0_wp 623 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 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(:,:) 624 632 DO jk = 2, jpk 625 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 626 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 627 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 633 DO jj = 1,jpj 634 DO ji = 1,jpi 635 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 636 ! 1 for jk = mikt 637 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 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) 642 END DO 643 END DO 628 644 END DO 629 ! Local depth and Inverse of the local depth of the water column at u- and v- points 630 ! ---------------------------------------------------------------------------------- 631 hu (:,:) = hu_a (:,:) 632 hv (:,:) = hv_a (:,:) 633 634 ! Inverse of the local depth 635 hur(:,:) = hur_a(:,:) 636 hvr(:,:) = hvr_a(:,:) 637 638 ! Local depth of the water column at t- points 639 ! -------------------------------------------- 640 ht(:,:) = 0. 641 DO jk = 1, jpkm1 642 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 645 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) 643 654 END DO 644 655 645 656 ! Write outputs 646 657 ! ============= 647 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 648 CALL iom_put( "cellthc" , fse3t_n (:,:,:) ) 649 CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 650 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 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(:,:,:) ) 663 IF( iom_use("e3tdef") ) & 664 CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 ) 651 665 652 666 ! write restart file 653 667 ! ================== 654 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 655 ! 656 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 657 ! 658 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 659 668 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 669 ! 670 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 671 ! 660 672 END SUBROUTINE dom_vvl_sf_swp 661 673 … … 671 683 !! - vertical interpolation: simple averaging 672 684 !!---------------------------------------------------------------------- 673 !! * Arguments 674 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 675 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 676 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors 677 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 678 !! * Local declarations 679 INTEGER :: ji, jj, jk ! dummy loop indices 680 LOGICAL :: l_is_orca ! local logical 681 !!---------------------------------------------------------------------- 682 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 683 ! 684 l_is_orca = .FALSE. 685 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations 686 687 SELECT CASE ( pout ) 688 ! ! ------------------------------------- ! 689 CASE( 'U' ) ! interpolation from T-point to U-point ! 690 ! ! ------------------------------------- ! 691 ! horizontal surface weighted interpolation 685 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated 686 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3 687 CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors 688 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 689 ! 690 INTEGER :: ji, jj, jk ! dummy loop indices 691 REAL(wp) :: zlnwd ! =1./0. when ln_wd = T/F 692 !!---------------------------------------------------------------------- 693 ! 694 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 695 ! 696 IF(ln_wd) THEN 697 zlnwd = 1.0_wp 698 ELSE 699 zlnwd = 0.0_wp 700 END IF 701 ! 702 SELECT CASE ( pout ) !== type of interpolation ==! 703 ! 704 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 692 705 DO jk = 1, jpk 693 706 DO jj = 1, jpjm1 694 707 DO ji = 1, fs_jpim1 ! vector opt. 695 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)&696 & * ( e1 2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &697 & + e1 2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )708 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 709 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 710 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 698 711 END DO 699 712 END DO 700 713 END DO 701 ! 702 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 703 ! boundary conditions 704 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. ) 714 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 705 715 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 706 ! ! ------------------------------------- ! 707 CASE( 'V' ) ! interpolation from T-point to V-point ! 708 ! ! ------------------------------------- ! 709 ! horizontal surface weighted interpolation 716 ! 717 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 710 718 DO jk = 1, jpk 711 719 DO jj = 1, jpjm1 712 720 DO ji = 1, fs_jpim1 ! vector opt. 713 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)&714 & * ( e1 2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &715 & + e1 2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )721 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 722 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 723 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 716 724 END DO 717 725 END DO 718 726 END DO 719 ! 720 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 721 ! boundary conditions 722 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. ) 727 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 723 728 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 724 ! ! ------------------------------------- ! 725 CASE( 'F' ) ! interpolation from U-point to F-point ! 726 ! ! ------------------------------------- ! 727 ! horizontal surface weighted interpolation 729 ! 730 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 728 731 DO jk = 1, jpk 729 732 DO jj = 1, jpjm1 730 733 DO ji = 1, fs_jpim1 ! vector opt. 731 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 732 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 733 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 734 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 735 & * r1_e1e2f(ji,jj) & 736 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 737 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 734 738 END DO 735 739 END DO 736 740 END DO 737 ! 738 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 739 ! boundary conditions 740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. ) 741 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 741 742 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 742 ! ! ------------------------------------- ! 743 CASE( 'W' ) ! interpolation from T-point to W-point ! 744 ! ! ------------------------------------- ! 745 ! vertical simple interpolation 743 ! 744 CASE( 'W' ) !* from T- to W-point : vertical simple mean 745 ! 746 746 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 747 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 747 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 748 !!gm BUG? use here wmask in case of ISF ? to be checked 748 749 DO jk = 2, jpk 749 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 750 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 751 END DO 752 ! ! -------------------------------------- ! 753 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 754 ! ! -------------------------------------- ! 755 ! vertical simple interpolation 750 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 751 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 752 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 753 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 754 END DO 755 ! 756 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean 757 ! 756 758 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 757 759 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 760 !!gm BUG? use here wumask in case of ISF ? to be checked 758 761 DO jk = 2, jpk 759 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 760 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 761 END DO 762 ! ! -------------------------------------- ! 763 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 764 ! ! -------------------------------------- ! 765 ! vertical simple interpolation 762 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 763 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 764 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 765 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 766 END DO 767 ! 768 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean 769 ! 766 770 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 767 771 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 772 !!gm BUG? use here wvmask in case of ISF ? to be checked 768 773 DO jk = 2, jpk 769 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 770 & + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 774 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 775 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 776 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 777 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 771 778 END DO 772 779 END SELECT 773 780 ! 774 775 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 776 781 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 782 ! 777 783 END SUBROUTINE dom_vvl_interpol 784 778 785 779 786 SUBROUTINE dom_vvl_rst( kt, cdrw ) … … 789 796 !! they are set to 0. 790 797 !!---------------------------------------------------------------------- 791 !! * Arguments792 798 INTEGER , INTENT(in) :: kt ! ocean time-step 793 799 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 794 ! ! * Local declarations795 INTEGER :: j k800 ! 801 INTEGER :: ji, jj, jk 796 802 INTEGER :: id1, id2, id3, id4, id5 ! local integers 797 803 !!---------------------------------------------------------------------- … … 804 810 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 805 811 ! 806 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )807 id2 = iom_varid( numror, ' fse3t_n', ldstop = .FALSE. )812 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 813 id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 808 814 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 809 815 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 810 id5 = iom_varid( numror, 'hdi f_lf', ldstop = .FALSE. )816 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 811 817 ! ! --------- ! 812 818 ! ! all cases ! 813 819 ! ! --------- ! 814 820 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 815 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 816 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 821 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 822 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 823 ! needed to restart if land processor not computed 824 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 825 WHERE ( tmask(:,:,:) == 0.0_wp ) 826 e3t_n(:,:,:) = e3t_0(:,:,:) 827 e3t_b(:,:,:) = e3t_0(:,:,:) 828 END WHERE 817 829 IF( neuler == 0 ) THEN 818 fse3t_b(:,:,:) = fse3t_n(:,:,:)830 e3t_b(:,:,:) = e3t_n(:,:,:) 819 831 ENDIF 820 832 ELSE IF( id1 > 0 ) THEN 821 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_bnot found in restart files'822 IF(lwp) write(numout,*) ' fse3t_b set equal to fse3t_n.'833 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 834 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 823 835 IF(lwp) write(numout,*) 'neuler is forced to 0' 824 fse3t_b(:,:,:) = fse3t_n(:,:,:) 836 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 837 e3t_n(:,:,:) = e3t_b(:,:,:) 838 neuler = 0 839 ELSE IF( id2 > 0 ) THEN 840 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 841 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 842 IF(lwp) write(numout,*) 'neuler is forced to 0' 843 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 844 e3t_b(:,:,:) = e3t_n(:,:,:) 825 845 neuler = 0 826 846 ELSE 827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'847 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 828 848 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 829 849 IF(lwp) write(numout,*) 'neuler is forced to 0' 830 DO jk =1,jpk831 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &832 & / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk)&833 & + e3t_0(:,:,jk)* (1._wp -tmask(:,:,jk))850 DO jk = 1, jpk 851 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 852 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 853 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 834 854 END DO 835 fse3t_b(:,:,:) = fse3t_n(:,:,:)855 e3t_b(:,:,:) = e3t_n(:,:,:) 836 856 neuler = 0 837 857 ENDIF … … 864 884 ! 865 885 ELSE !* Initialize at "rest" 866 fse3t_b(:,:,:) = e3t_0(:,:,:)867 fse3t_n(:,:,:) = e3t_0(:,:,:)886 e3t_b(:,:,:) = e3t_0(:,:,:) 887 e3t_n(:,:,:) = e3t_0(:,:,:) 868 888 sshn(:,:) = 0.0_wp 889 890 IF( ln_wd ) THEN 891 DO jj = 1, jpj 892 DO ji = 1, jpi 893 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 894 e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 895 e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 896 e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 897 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 898 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 899 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 900 ENDIF 901 ENDDO 902 ENDDO 903 END IF 904 869 905 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 870 906 tilde_e3t_b(:,:,:) = 0.0_wp … … 873 909 END IF 874 910 ENDIF 875 911 ! 876 912 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 877 913 ! ! =================== … … 880 916 ! ! all cases ! 881 917 ! ! --------- ! 882 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )883 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_n', fse3t_n(:,:,:) )918 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) ) 919 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) ) 884 920 ! ! ----------------------- ! 885 921 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 893 929 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 894 930 ENDIF 895 896 ENDIF 931 ! 932 ENDIF 933 ! 897 934 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 898 935 ! 899 936 END SUBROUTINE dom_vvl_rst 900 937 … … 907 944 !! for vertical coordinate 908 945 !!---------------------------------------------------------------------- 909 INTEGER :: ioptio 910 INTEGER :: ios 911 946 INTEGER :: ioptio, ios 947 !! 912 948 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 913 &ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &914 &rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe949 & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & 950 & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 915 951 !!---------------------------------------------------------------------- 916 952 ! 917 953 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist : 918 954 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 919 955 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 920 956 ! 921 957 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run 922 958 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 923 959 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) 924 960 IF(lwm) WRITE ( numond, nam_vvl ) 925 961 ! 926 962 IF(lwp) THEN ! Namelist print 927 963 WRITE(numout,*) … … 956 992 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 957 993 ENDIF 958 994 ! 959 995 ioptio = 0 ! Parameter control 960 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.961 IF( ln_vvl_zstar ) 962 IF( ln_vvl_ztilde ) 963 IF( ln_vvl_layer ) 964 996 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. 997 IF( ln_vvl_zstar ) ioptio = ioptio + 1 998 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 999 IF( ln_vvl_layer ) ioptio = ioptio + 1 1000 ! 965 1001 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 966 1002 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 1003 ! 967 1004 IF(lwp) THEN ! Print the choice 968 1005 WRITE(numout,*) … … 975 1012 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 976 1013 ENDIF 977 1014 ! 978 1015 #if defined key_agrif 979 IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )1016 IF(.NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' ) 980 1017 #endif 981 1018 ! 982 1019 END SUBROUTINE dom_vvl_ctl 983 984 SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )985 !!---------------------------------------------------------------------986 !! *** ROUTINE dom_vvl_orca_fix ***987 !!988 !! ** Purpose : Correct surface weighted, horizontally interpolated,989 !! scale factors at locations that have been individually990 !! modified in domhgr. Such modifications break the991 !! relationship between e12t and e1u*e2u etc.992 !! Recompute some scale factors ignoring the modified metric.993 !!----------------------------------------------------------------------994 !! * Arguments995 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated996 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3997 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors998 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'999 !! * Local declarations1000 INTEGER :: ji, jj, jk ! dummy loop indices1001 INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices1002 !! acc1003 !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for1004 !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations1005 !!1006 ! ! =====================1007 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration1008 ! ! =====================1009 !! acc1010 IF( nn_cla == 0 ) THEN1011 !1012 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)1013 ij0 = 102 ; ij1 = 1021014 DO jk = 1, jpkm11015 DO jj = mj0(ij0), mj1(ij1)1016 DO ji = mi0(ii0), mi1(ii1)1017 SELECT CASE ( pout )1018 CASE( 'U' )1019 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1020 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1021 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1022 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1023 CASE( 'F' )1024 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1025 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1026 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1027 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1028 END SELECT1029 END DO1030 END DO1031 END DO1032 !1033 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)1034 ij0 = 88 ; ij1 = 881035 DO jk = 1, jpkm11036 DO jj = mj0(ij0), mj1(ij1)1037 DO ji = mi0(ii0), mi1(ii1)1038 SELECT CASE ( pout )1039 CASE( 'U' )1040 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1041 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1042 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1043 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1044 CASE( 'V' )1045 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1046 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1047 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1048 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1049 CASE( 'F' )1050 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1051 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1052 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1053 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1054 END SELECT1055 END DO1056 END DO1057 END DO1058 ENDIF1059 1060 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)1061 ij0 = 116 ; ij1 = 1161062 DO jk = 1, jpkm11063 DO jj = mj0(ij0), mj1(ij1)1064 DO ji = mi0(ii0), mi1(ii1)1065 SELECT CASE ( pout )1066 CASE( 'U' )1067 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1068 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1069 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1070 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1071 CASE( 'F' )1072 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1073 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1074 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1075 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1076 END SELECT1077 END DO1078 END DO1079 END DO1080 ENDIF1081 !1082 ! ! =====================1083 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration1084 ! ! =====================1085 !1086 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified)1087 ij0 = 200 ; ij1 = 2001088 DO jk = 1, jpkm11089 DO jj = mj0(ij0), mj1(ij1)1090 DO ji = mi0(ii0), mi1(ii1)1091 SELECT CASE ( pout )1092 CASE( 'U' )1093 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1094 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1095 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1096 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1097 CASE( 'F' )1098 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1099 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1100 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1101 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1102 END SELECT1103 END DO1104 END DO1105 END DO1106 !1107 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)1108 ij0 = 208 ; ij1 = 2081109 DO jk = 1, jpkm11110 DO jj = mj0(ij0), mj1(ij1)1111 DO ji = mi0(ii0), mi1(ii1)1112 SELECT CASE ( pout )1113 CASE( 'U' )1114 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1115 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1116 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1117 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1118 CASE( 'F' )1119 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1120 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1121 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1122 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1123 END SELECT1124 END DO1125 END DO1126 END DO1127 !1128 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)1129 ij0 = 124 ; ij1 = 1251130 DO jk = 1, jpkm11131 DO jj = mj0(ij0), mj1(ij1)1132 DO ji = mi0(ii0), mi1(ii1)1133 SELECT CASE ( pout )1134 CASE( 'V' )1135 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1136 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1137 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1138 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1139 END SELECT1140 END DO1141 END DO1142 END DO1143 !1144 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]1145 ij0 = 124 ; ij1 = 1251146 DO jk = 1, jpkm11147 DO jj = mj0(ij0), mj1(ij1)1148 DO ji = mi0(ii0), mi1(ii1)1149 SELECT CASE ( pout )1150 CASE( 'V' )1151 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1152 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1153 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1154 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1155 END SELECT1156 END DO1157 END DO1158 END DO1159 !1160 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)1161 ij0 = 124 ; ij1 = 1251162 DO jk = 1, jpkm11163 DO jj = mj0(ij0), mj1(ij1)1164 DO ji = mi0(ii0), mi1(ii1)1165 SELECT CASE ( pout )1166 CASE( 'V' )1167 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1168 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1169 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1170 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1171 END SELECT1172 END DO1173 END DO1174 END DO1175 !1176 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)1177 ij0 = 124 ; ij1 = 1251178 DO jk = 1, jpkm11179 DO jj = mj0(ij0), mj1(ij1)1180 DO ji = mi0(ii0), mi1(ii1)1181 SELECT CASE ( pout )1182 CASE( 'V' )1183 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1184 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1185 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1186 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1187 END SELECT1188 END DO1189 END DO1190 END DO1191 !1192 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)1193 ij0 = 141 ; ij1 = 1421194 DO jk = 1, jpkm11195 DO jj = mj0(ij0), mj1(ij1)1196 DO ji = mi0(ii0), mi1(ii1)1197 SELECT CASE ( pout )1198 CASE( 'V' )1199 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1200 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1201 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1202 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1203 END SELECT1204 END DO1205 END DO1206 END DO1207 !1208 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)1209 ij0 = 141 ; ij1 = 1421210 DO jk = 1, jpkm11211 DO jj = mj0(ij0), mj1(ij1)1212 DO ji = mi0(ii0), mi1(ii1)1213 SELECT CASE ( pout )1214 CASE( 'V' )1215 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1216 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1217 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1218 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1219 END SELECT1220 END DO1221 END DO1222 END DO1223 ENDIF1224 ! ! =====================1225 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration1226 ! ! =====================1227 !1228 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)1229 ij0 = 327 ; ij1 = 3271230 DO jk = 1, jpkm11231 DO jj = mj0(ij0), mj1(ij1)1232 DO ji = mi0(ii0), mi1(ii1)1233 SELECT CASE ( pout )1234 CASE( 'U' )1235 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1236 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1237 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1238 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1239 CASE( 'F' )1240 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1241 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1242 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1243 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1244 END SELECT1245 END DO1246 END DO1247 END DO1248 !1249 ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)1250 ij0 = 343 ; ij1 = 3431251 DO jk = 1, jpkm11252 DO jj = mj0(ij0), mj1(ij1)1253 DO ji = mi0(ii0), mi1(ii1)1254 SELECT CASE ( pout )1255 CASE( 'U' )1256 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1257 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1258 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1259 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1260 CASE( 'F' )1261 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1262 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1263 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1264 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1265 END SELECT1266 END DO1267 END DO1268 END DO1269 !1270 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)1271 ij0 = 232 ; ij1 = 2321272 DO jk = 1, jpkm11273 DO jj = mj0(ij0), mj1(ij1)1274 DO ji = mi0(ii0), mi1(ii1)1275 SELECT CASE ( pout )1276 CASE( 'U' )1277 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1278 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1279 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1280 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1281 CASE( 'F' )1282 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1283 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1284 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1285 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1286 END SELECT1287 END DO1288 END DO1289 END DO1290 !1291 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)1292 ij0 = 232 ; ij1 = 2321293 DO jk = 1, jpkm11294 DO jj = mj0(ij0), mj1(ij1)1295 DO ji = mi0(ii0), mi1(ii1)1296 SELECT CASE ( pout )1297 CASE( 'U' )1298 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1299 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1300 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1301 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1302 CASE( 'F' )1303 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1304 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1305 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1306 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1307 END SELECT1308 END DO1309 END DO1310 END DO1311 !1312 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)1313 ij0 = 270 ; ij1 = 2701314 DO jk = 1, jpkm11315 DO jj = mj0(ij0), mj1(ij1)1316 DO ji = mi0(ii0), mi1(ii1)1317 SELECT CASE ( pout )1318 CASE( 'U' )1319 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1320 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1321 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1322 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1323 CASE( 'F' )1324 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1325 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1326 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1327 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1328 END SELECT1329 END DO1330 END DO1331 END DO1332 !1333 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)1334 ij0 = 232 ; ij1 = 2331335 DO jk = 1, jpkm11336 DO jj = mj0(ij0), mj1(ij1)1337 DO ji = mi0(ii0), mi1(ii1)1338 SELECT CASE ( pout )1339 CASE( 'V' )1340 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1341 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1342 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1343 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1344 END SELECT1345 END DO1346 END DO1347 END DO1348 !1349 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)1350 ij0 = 276 ; ij1 = 2761351 DO jk = 1, jpkm11352 DO jj = mj0(ij0), mj1(ij1)1353 DO ji = mi0(ii0), mi1(ii1)1354 SELECT CASE ( pout )1355 CASE( 'V' )1356 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1357 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1358 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1359 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1360 END SELECT1361 END DO1362 END DO1363 END DO1364 ENDIF1365 END SUBROUTINE dom_vvl_orca_fix1366 1020 1367 1021 !!====================================================================== 1368 1022 END MODULE domvvl 1369 1370 1371
Note: See TracChangeset
for help on using the changeset viewer.
