- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5506 r6808 10 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 11 11 !!---------------------------------------------------------------------- 12 !! 'key_vvl' variable volume 13 !!---------------------------------------------------------------------- 12 14 13 !!---------------------------------------------------------------------- 15 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness … … 19 18 !! dom_vvl_rst : read/write restart file 20 19 !! dom_vvl_ctl : Check the vvl options 21 !! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors22 !! : to account for manual changes to e[1,2][u,v] in some Straits23 20 !!---------------------------------------------------------------------- 24 !! * Modules used25 21 USE oce ! ocean dynamics and tracers 22 USE phycst ! physical constant 26 23 USE dom_oce ! ocean space and time domain 27 24 USE sbc_oce ! ocean surface boundary condition 25 USE wet_dry ! wetting and drying 26 USE restart ! ocean restart 27 ! 28 28 USE in_out_manager ! I/O manager 29 29 USE iom ! I/O manager library 30 USE restart ! ocean restart31 30 USE lib_mpp ! distributed memory computing library 32 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 37 36 PRIVATE 38 37 39 !! * Routine accessibility40 38 PUBLIC dom_vvl_init ! called by domain.F90 41 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 42 40 PUBLIC dom_vvl_sf_swp ! called by step.F90 43 41 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 44 PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol 45 46 !!* Namelist nam_vvl 47 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 53 ! ! conservation: not used yet 54 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 55 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 56 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 57 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 58 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 59 60 !! * Module variables 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 62 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 66 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 67 63 68 64 !! * Substitutions 69 # include "domzgr_substitute.h90"70 65 # include "vectopt_loop_substitute.h90" 71 66 !!---------------------------------------------------------------------- 72 !! NEMO/OPA 3. 3 , NEMO-Consortium (2010)67 !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 73 68 !! $Id$ 74 69 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 75 70 !!---------------------------------------------------------------------- 76 77 71 CONTAINS 78 72 … … 81 75 !! *** FUNCTION dom_vvl_alloc *** 82 76 !!---------------------------------------------------------------------- 83 IF( ln_vvl_zstar ) dom_vvl_alloc = 077 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 84 78 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 85 79 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & … … 88 82 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 89 83 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 90 un_td = 0. 0_wp91 vn_td = 0. 0_wp84 un_td = 0._wp 85 vn_td = 0._wp 92 86 ENDIF 93 87 IF( ln_vvl_ztilde ) THEN … … 96 90 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 97 91 ENDIF 98 92 ! 99 93 END FUNCTION dom_vvl_alloc 100 94 … … 110 104 !! - interpolate scale factors 111 105 !! 112 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b)113 !! - Regrid: fse3(u/v)_n114 !! fse3(u/v)_b115 !! fse3w_n116 !! fse3(u/v)w_b117 !! fse3(u/v)w_n118 !! 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 119 113 !! - h(t/u/v)_0 120 114 !! - frq_rst_e3t and frq_rst_hdv … … 122 116 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 123 117 !!---------------------------------------------------------------------- 124 USE phycst, ONLY : rpi, rsmall, rad 125 !! * Local declarations 126 INTEGER :: ji,jj,jk 118 INTEGER :: ji, jj, jk 127 119 INTEGER :: ii0, ii1, ij0, ij1 128 120 REAL(wp):: zcoef 129 121 !!---------------------------------------------------------------------- 130 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 131 122 ! 123 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 124 ! 132 125 IF(lwp) WRITE(numout,*) 133 126 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 134 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 135 136 ! choose vertical coordinate (z_star, z_tilde or layer) 137 ! ========================== 138 CALL dom_vvl_ctl 139 140 ! Allocate module arrays 141 ! ====================== 128 ! 129 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 130 ! 131 ! ! Allocate module arrays 142 132 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 143 144 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 145 ! ============================================================================ 133 ! 134 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 146 135 CALL dom_vvl_rst( nit000, 'READ' ) 147 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 148 149 ! Reconstruction of all vertical scale factors at now and before time steps 150 ! ============================================================================= 151 ! Horizontal scale factor interpolations 152 ! -------------------------------------- 153 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 154 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 155 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 156 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 157 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 158 ! Vertical scale factor interpolations 159 ! ------------------------------------ 160 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 161 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 162 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 163 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 164 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 165 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 166 ! t- and w- points depth 167 ! ---------------------- 168 ! set the isf depth as it is in the initial step 169 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 170 fsdepw_n(:,:,1) = 0.0_wp 171 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 172 fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 173 fsdepw_b(:,:,1) = 0.0_wp 174 175 DO jk = 2, jpk 136 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 137 ! 138 ! !== Set of all other vertical scale factors ==! (now and before) 139 ! ! Horizontal interpolation of e3t 140 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 141 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 142 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 143 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 144 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 145 ! ! Vertical interpolation of e3t,u,v 146 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 147 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 148 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 149 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 150 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 151 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 152 ! 153 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 154 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 155 gdepw_n(:,:,1) = 0.0_wp 156 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 157 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 158 gdepw_b(:,:,1) = 0.0_wp 159 DO jk = 2, jpk ! vertical sum 176 160 DO jj = 1,jpj 177 161 DO ji = 1,jpi 178 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 179 ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 180 ! 0.5 where jk = mikt 181 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 182 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 183 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 184 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 185 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 186 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 187 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 188 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 162 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 163 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 164 ! ! 0.5 where jk = mikt 165 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 166 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 167 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 168 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 169 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 170 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 171 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 172 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 173 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 189 174 END DO 190 175 END DO 191 176 END DO 192 193 ! Before depth and Inverse of the local depth of the water column at u- and v- points 194 ! ----------------------------------------------------------------------------------- 195 hu_b(:,:) = 0. 196 hv_b(:,:) = 0. 197 DO jk = 1, jpkm1 198 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 199 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) 200 190 END DO 201 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 202 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 203 204 ! Restoring frequencies for z_tilde coordinate 205 ! ============================================ 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) 206 199 IF( ln_vvl_ztilde ) THEN 207 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 208 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 209 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 210 IF( ln_vvl_ztilde_as_zstar ) THEN 211 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 212 frq_rst_e3t(:,:) = 0.0_wp 213 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 214 209 ENDIF 215 IF ( ln_vvl_zstar_at_eqtor ) THEN 210 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 216 211 DO jj = 1, jpj 217 212 DO ji = 1, jpi 213 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 218 214 IF( ABS(gphit(ji,jj)) >= 6.) THEN 219 215 ! values outside the equatorial band and transition zone (ztilde) 220 216 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 221 217 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 222 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 218 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 223 219 ! values inside the equatorial band (ztilde as zstar) 224 220 frq_rst_e3t(ji,jj) = 0.0_wp 225 221 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 226 ELSE 227 ! 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) 228 224 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 229 225 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & … … 236 232 END DO 237 233 END DO 238 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 239 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 240 236 ij0 = 128 ; ij1 = 135 ; 241 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp … … 244 240 ENDIF 245 241 ENDIF 246 242 ! 247 243 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 248 244 ! 249 245 END SUBROUTINE dom_vvl_init 250 246 … … 268 264 !! - tilde_e3t_a: after increment of vertical scale factor 269 265 !! in z_tilde case 270 !! - fse3(t/u/v)_a266 !! - e3(t/u/v)_a 271 267 !! 272 268 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 273 269 !!---------------------------------------------------------------------- 274 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 275 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 276 !! * Arguments 277 INTEGER, INTENT( in ) :: kt ! time step 278 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 279 !! * Local declarations 280 INTEGER :: ji, jj, jk ! dummy loop indices 281 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 282 REAL(wp) :: z2dt ! temporary scalars 283 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 284 LOGICAL :: ll_do_bclinic ! temporary logical 285 !!---------------------------------------------------------------------- 286 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 287 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 288 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 289 290 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 291 289 IF(lwp) WRITE(numout,*) 292 290 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' … … 296 294 ll_do_bclinic = .TRUE. 297 295 IF( PRESENT(kcall) ) THEN 298 IF ( kcall == 2 .AND. ln_vvl_ztilde )ll_do_bclinic = .FALSE.296 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 299 297 ENDIF 300 298 … … 302 300 ! After acale factors at t-points ! 303 301 ! ******************************* ! 304 305 302 ! ! --------------------------------------------- ! 306 303 ! ! z_star coordinate and barotropic z-tilde part ! 307 304 ! ! --------------------------------------------- ! 308 305 ! 309 306 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 310 307 DO jk = 1, jpkm1 311 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)312 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) 313 310 END DO 314 311 ! 315 312 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 316 313 ! ! ------baroclinic part------ ! 317 318 314 ! I - initialization 319 315 ! ================== … … 321 317 ! 1 - barotropic divergence 322 318 ! ------------------------- 323 zhdiv(:,:) = 0. 324 zht(:,:) = 0. 319 zhdiv(:,:) = 0._wp 320 zht(:,:) = 0._wp 325 321 DO jk = 1, jpkm1 326 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)327 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)322 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 323 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 328 324 END DO 329 325 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 332 328 ! -------------------------------------------------- 333 329 IF( ln_vvl_ztilde ) THEN 334 IF( kt .GT.nit000 ) THEN330 IF( kt > nit000 ) THEN 335 331 DO jk = 1, jpkm1 336 332 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 337 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )333 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 338 334 END DO 339 335 ENDIF 340 END 336 ENDIF 341 337 342 338 ! II - after z_tilde increments of vertical scale factors 343 339 ! ======================================================= 344 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 345 341 346 342 ! 1 - High frequency divergence term … … 348 344 IF( ln_vvl_ztilde ) THEN ! z_tilde case 349 345 DO jk = 1, jpkm1 350 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) ) 351 347 END DO 352 348 ELSE ! layer case 353 349 DO jk = 1, jpkm1 354 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)355 END DO 356 END 350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 351 END DO 352 ENDIF 357 353 358 354 ! 2 - Restoring term (z-tilde case only) … … 362 358 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 363 359 END DO 364 END 360 ENDIF 365 361 366 362 ! 3 - Thickness diffusion term 367 363 ! ---------------------------- 368 zwu(:,:) = 0.0_wp 369 zwv(:,:) = 0.0_wp 370 ! a - first derivative: diffusive fluxes 371 DO jk = 1, jpkm1 364 zwu(:,:) = 0._wp 365 zwv(:,:) = 0._wp 366 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 372 367 DO jj = 1, jpjm1 373 368 DO ji = 1, fs_jpim1 ! vector opt. 374 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj)&375 &* ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) )376 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj)&377 &* ( 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) ) 378 373 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 379 374 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) … … 381 376 END DO 382 377 END DO 383 ! b - correction for last oceanic u-v points 384 DO jj = 1, jpj 378 DO jj = 1, jpj ! b - correction for last oceanic u-v points 385 379 DO ji = 1, jpi 386 380 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) … … 388 382 END DO 389 383 END DO 390 ! c - second derivative: divergence of diffusive fluxes 391 DO jk = 1, jpkm1 384 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 392 385 DO jj = 2, jpjm1 393 386 DO ji = fs_2, fs_jpim1 ! vector opt. 394 387 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 395 388 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 396 & ) * r1_e1 2t(ji,jj)389 & ) * r1_e1e2t(ji,jj) 397 390 END DO 398 391 END DO 399 392 END DO 400 ! d - thickness diffusion transport: boundary conditions401 ! (stored for tracer advction and continuity equation)393 ! ! d - thickness diffusion transport: boundary conditions 394 ! (stored for tracer advction and continuity equation) 402 395 CALL lbc_lnk( un_td , 'U' , -1._wp) 403 396 CALL lbc_lnk( vn_td , 'V' , -1._wp) … … 417 410 ! Maximum deformation control 418 411 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 419 ze3t(:,:,jpk) = 0. 0_wp412 ze3t(:,:,jpk) = 0._wp 420 413 DO jk = 1, jpkm1 421 414 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) … … 426 419 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 427 420 ! - ML - test: for the moment, stop simulation for too large e3_t variations 428 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 429 422 IF( lk_mpp ) THEN 430 423 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) … … 469 462 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 470 463 DO jk = 1, jpkm1 471 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) 472 465 END DO 473 466 … … 477 470 ! ! ---baroclinic part--------- ! 478 471 DO jk = 1, jpkm1 479 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)472 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 480 473 END DO 481 474 ENDIF … … 492 485 zht(:,:) = 0.0_wp 493 486 DO jk = 1, jpkm1 494 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)487 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 495 488 END DO 496 489 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 497 490 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 498 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 499 492 ! 500 493 zht(:,:) = 0.0_wp 501 494 DO jk = 1, jpkm1 502 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)495 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 503 496 END DO 504 497 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 505 498 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 506 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 507 500 ! 508 501 zht(:,:) = 0.0_wp 509 502 DO jk = 1, jpkm1 510 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)503 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 511 504 END DO 512 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 513 506 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 514 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 515 508 ! 516 509 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) … … 531 524 ! *********************************** ! 532 525 533 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )534 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' ) 535 528 536 529 ! *********************************** ! … … 538 531 ! *********************************** ! 539 532 540 hu_a(:,:) = 0._wp ! Ocean depth at U-points541 hv_a(:,:) = 0._wp ! Ocean depth at V-points542 DO jk = 1, jpkm1543 hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)544 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) 545 538 END DO 546 539 ! ! Inverse of the local depth 547 hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 548 hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 549 550 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 551 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 552 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 ! 553 547 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') 554 548 ! 555 549 END SUBROUTINE dom_vvl_sf_nxt 556 550 … … 568 562 !! - recompute depths and water height fields 569 563 !! 570 !! ** 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 571 565 !! - Recompute: 572 !! fse3(u/v)_b573 !! fse3w_n574 !! fse3(u/v)w_b575 !! fse3(u/v)w_n576 !! 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 577 571 !! h(u/v) and h(u/v)r 578 572 !! … … 580 574 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 581 575 !!---------------------------------------------------------------------- 582 !! * Arguments 583 INTEGER, INTENT( in ) :: kt ! time step 584 !! * Local declarations 585 INTEGER :: ji,jj,jk ! dummy loop indices 586 REAL(wp) :: zcoef 587 !!---------------------------------------------------------------------- 588 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 ! 589 584 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 590 585 ! … … 597 592 ! Time filter and swap of scale factors 598 593 ! ===================================== 599 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.594 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 600 595 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 601 596 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 607 602 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 608 603 ENDIF 609 fsdept_b(:,:,:) = fsdept_n(:,:,:)610 fsdepw_b(:,:,:) = fsdepw_n(:,:,:)611 612 fse3t_n(:,:,:) = fse3t_a(:,:,:)613 fse3u_n(:,:,:) = fse3u_a(:,:,:)614 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(:,:,:) 615 610 616 611 ! Compute all missing vertical scale factor and depths … … 618 613 ! Horizontal scale factor interpolations 619 614 ! -------------------------------------- 620 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt615 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 621 616 ! - JC - hu_b, hv_b, hur_b, hvr_b also 622 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 617 618 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 619 623 620 ! Vertical scale factor interpolations 624 ! ------------------------------------ 625 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 626 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 627 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 628 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W' ) 629 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 630 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 631 ! t- and w- points depth 632 ! ---------------------- 633 ! set the isf depth as it is in the initial step 634 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 635 fsdepw_n(:,:,1) = 0.0_wp 636 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 637 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(:,:) 638 632 DO jk = 2, jpk 639 633 DO jj = 1,jpj … … 642 636 ! 1 for jk = mikt 643 637 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 644 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)645 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) &646 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))647 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)638 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 639 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 640 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 641 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 648 642 END DO 649 643 END DO 650 644 END DO 651 645 652 ! Local depth and Inverse of the local depth of the water column at u- and v- points 653 ! ---------------------------------------------------------------------------------- 654 hu (:,:) = hu_a (:,:) 655 hv (:,:) = hv_a (:,:) 656 657 ! Inverse of the local depth 658 hur(:,:) = hur_a(:,:) 659 hvr(:,:) = hvr_a(:,:) 660 661 ! Local depth of the water column at t- points 662 ! -------------------------------------------- 663 ht(:,:) = 0. 664 DO jk = 1, jpkm1 665 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 646 ! Local depth and Inverse of the local depth of the water 647 ! ------------------------------------------------------- 648 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 649 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 650 ! 651 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 652 DO jk = 2, jpkm1 653 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 666 654 END DO 667 655 668 656 ! Write outputs 669 657 ! ============= 670 CALL iom_put( "e3t" , fse3t_n(:,:,:) )671 CALL iom_put( "e3u" , fse3u_n(:,:,:) )672 CALL iom_put( "e3v" , fse3v_n(:,:,:) )673 CALL iom_put( "e3w" , fse3w_n(:,:,:) )674 CALL iom_put( "tpt_dep" , fsde3w_n(:,:,:) )658 CALL iom_put( "e3t", e3t_n(:,:,:) ) 659 CALL iom_put( "e3u", e3u_n(:,:,:) ) 660 CALL iom_put( "e3v", e3v_n(:,:,:) ) 661 CALL iom_put( "e3w", e3w_n(:,:,:) ) 662 CALL iom_put( "tpt_dep", gde3w_n(:,:,:) ) 675 663 IF( iom_use("e3tdef") ) & 676 CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100* tmask(:,:,:) ) ** 2 )664 CALL iom_put( "e3tdef", ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100. * tmask(:,:,:) ) ** 2 ) 677 665 678 666 ! write restart file 679 667 ! ================== 680 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )681 ! 682 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp')683 668 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 669 ! 670 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 671 ! 684 672 END SUBROUTINE dom_vvl_sf_swp 685 673 … … 695 683 !! - vertical interpolation: simple averaging 696 684 !!---------------------------------------------------------------------- 697 !! * Arguments 698 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 699 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 700 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors 701 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 702 !! * Local declarations 703 INTEGER :: ji, jj, jk ! dummy loop indices 704 LOGICAL :: l_is_orca ! local logical 705 !!---------------------------------------------------------------------- 706 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 707 ! 708 l_is_orca = .FALSE. 709 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations 710 711 SELECT CASE ( pout ) 712 ! ! ------------------------------------- ! 713 CASE( 'U' ) ! interpolation from T-point to U-point ! 714 ! ! ------------------------------------- ! 715 ! 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 716 705 DO jk = 1, jpk 717 706 DO jj = 1, jpjm1 718 707 DO ji = 1, fs_jpim1 ! vector opt. 719 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)&720 & * ( e1 2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &721 & + 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) ) ) 722 711 END DO 723 712 END DO 724 713 END DO 725 !726 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )727 ! boundary conditions728 714 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 729 715 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 730 ! ! ------------------------------------- ! 731 CASE( 'V' ) ! interpolation from T-point to V-point ! 732 ! ! ------------------------------------- ! 733 ! horizontal surface weighted interpolation 716 ! 717 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 734 718 DO jk = 1, jpk 735 719 DO jj = 1, jpjm1 736 720 DO ji = 1, fs_jpim1 ! vector opt. 737 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)&738 & * ( e1 2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &739 & + 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) ) ) 740 724 END DO 741 725 END DO 742 726 END DO 743 !744 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )745 ! boundary conditions746 727 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 747 728 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 748 ! ! ------------------------------------- ! 749 CASE( 'F' ) ! interpolation from U-point to F-point ! 750 ! ! ------------------------------------- ! 751 ! horizontal surface weighted interpolation 729 ! 730 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 752 731 DO jk = 1, jpk 753 732 DO jj = 1, jpjm1 754 733 DO ji = 1, fs_jpim1 ! vector opt. 755 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 756 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 757 & + 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) ) ) 758 738 END DO 759 739 END DO 760 740 END DO 761 !762 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )763 ! boundary conditions764 741 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 765 742 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 766 ! ! ------------------------------------- ! 767 CASE( 'W' ) ! interpolation from T-point to W-point ! 768 ! ! ------------------------------------- ! 769 ! vertical simple interpolation 743 ! 744 CASE( 'W' ) !* from T- to W-point : vertical simple mean 745 ! 770 746 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 771 ! - 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 772 749 DO jk = 2, jpk 773 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 774 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 775 END DO 776 ! ! -------------------------------------- ! 777 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 778 ! ! -------------------------------------- ! 779 ! 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 ! 780 758 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 781 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 782 761 DO jk = 2, jpk 783 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 784 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 785 END DO 786 ! ! -------------------------------------- ! 787 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 788 ! ! -------------------------------------- ! 789 ! 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 ! 790 770 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 791 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 792 773 DO jk = 2, jpk 793 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 794 & + 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 ) ) 795 778 END DO 796 779 END SELECT 797 780 ! 798 799 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 800 781 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 782 ! 801 783 END SUBROUTINE dom_vvl_interpol 784 802 785 803 786 SUBROUTINE dom_vvl_rst( kt, cdrw ) … … 813 796 !! they are set to 0. 814 797 !!---------------------------------------------------------------------- 815 !! * Arguments816 798 INTEGER , INTENT(in) :: kt ! ocean time-step 817 799 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 818 ! ! * Local declarations819 INTEGER :: j k800 ! 801 INTEGER :: ji, jj, jk 820 802 INTEGER :: id1, id2, id3, id4, id5 ! local integers 821 803 !!---------------------------------------------------------------------- … … 828 810 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 829 811 ! 830 id1 = iom_varid( numror, ' fse3t_b', ldstop = .FALSE. )831 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. ) 832 814 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 833 815 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) … … 837 819 ! ! --------- ! 838 820 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 839 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )840 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(:,:,:) ) 841 823 ! needed to restart if land processor not computed 842 IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'824 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 843 825 WHERE ( tmask(:,:,:) == 0.0_wp ) 844 fse3t_n(:,:,:) = e3t_0(:,:,:)845 fse3t_b(:,:,:) = e3t_0(:,:,:)826 e3t_n(:,:,:) = e3t_0(:,:,:) 827 e3t_b(:,:,:) = e3t_0(:,:,:) 846 828 END WHERE 847 829 IF( neuler == 0 ) THEN 848 fse3t_b(:,:,:) = fse3t_n(:,:,:)830 e3t_b(:,:,:) = e3t_n(:,:,:) 849 831 ENDIF 850 832 ELSE IF( id1 > 0 ) THEN 851 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'852 IF(lwp) write(numout,*) ' fse3t_n set equal to fse3t_b.'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.' 853 835 IF(lwp) write(numout,*) 'neuler is forced to 0' 854 CALL iom_get( numror, jpdom_autoglo, ' fse3t_b', fse3t_b(:,:,:) )855 fse3t_n(:,:,:) = fse3t_b(:,:,:)836 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) ) 837 e3t_n(:,:,:) = e3t_b(:,:,:) 856 838 neuler = 0 857 839 ELSE IF( id2 > 0 ) THEN 858 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'859 IF(lwp) write(numout,*) ' fse3t_b set equal to fse3t_n.'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.' 860 842 IF(lwp) write(numout,*) 'neuler is forced to 0' 861 CALL iom_get( numror, jpdom_autoglo, ' fse3t_n', fse3t_n(:,:,:) )862 fse3t_b(:,:,:) = fse3t_n(:,:,:)843 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) ) 844 e3t_b(:,:,:) = e3t_n(:,:,:) 863 845 neuler = 0 864 846 ELSE 865 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' 866 848 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 867 849 IF(lwp) write(numout,*) 'neuler is forced to 0' 868 DO jk =1,jpk869 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &870 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)&871 & + 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)) 872 854 END DO 873 fse3t_b(:,:,:) = fse3t_n(:,:,:)855 e3t_b(:,:,:) = e3t_n(:,:,:) 874 856 neuler = 0 875 857 ENDIF … … 902 884 ! 903 885 ELSE !* Initialize at "rest" 904 fse3t_b(:,:,:) = e3t_0(:,:,:)905 fse3t_n(:,:,:) = e3t_0(:,:,:)886 e3t_b(:,:,:) = e3t_0(:,:,:) 887 e3t_n(:,:,:) = e3t_0(:,:,:) 906 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 907 905 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 908 906 tilde_e3t_b(:,:,:) = 0.0_wp … … 911 909 END IF 912 910 ENDIF 913 911 ! 914 912 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 915 913 ! ! =================== … … 918 916 ! ! all cases ! 919 917 ! ! --------- ! 920 CALL iom_rstput( kt, nitrst, numrow, ' fse3t_b', fse3t_b(:,:,:) )921 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(:,:,:) ) 922 920 ! ! ----------------------- ! 923 921 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 931 929 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 932 930 ENDIF 933 934 ENDIF 931 ! 932 ENDIF 933 ! 935 934 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 936 935 ! 937 936 END SUBROUTINE dom_vvl_rst 938 937 … … 945 944 !! for vertical coordinate 946 945 !!---------------------------------------------------------------------- 947 INTEGER :: ioptio 948 INTEGER :: ios 949 946 INTEGER :: ioptio, ios 947 !! 950 948 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 951 &ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , &952 &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 953 951 !!---------------------------------------------------------------------- 954 952 ! 955 953 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist : 956 954 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 957 955 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 958 956 ! 959 957 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run 960 958 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 961 959 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) 962 960 IF(lwm) WRITE ( numond, nam_vvl ) 963 961 ! 964 962 IF(lwp) THEN ! Namelist print 965 963 WRITE(numout,*) … … 994 992 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 995 993 ENDIF 996 994 ! 997 995 ioptio = 0 ! Parameter control 998 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.999 IF( ln_vvl_zstar ) 1000 IF( ln_vvl_ztilde ) 1001 IF( ln_vvl_layer ) 1002 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 ! 1003 1001 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1004 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1005 1002 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 1003 ! 1006 1004 IF(lwp) THEN ! Print the choice 1007 1005 WRITE(numout,*) … … 1014 1012 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 1015 1013 ENDIF 1016 1014 ! 1017 1015 #if defined key_agrif 1018 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' ) 1019 1017 #endif 1020 1018 ! 1021 1019 END SUBROUTINE dom_vvl_ctl 1022 1023 SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )1024 !!---------------------------------------------------------------------1025 !! *** ROUTINE dom_vvl_orca_fix ***1026 !!1027 !! ** Purpose : Correct surface weighted, horizontally interpolated,1028 !! scale factors at locations that have been individually1029 !! modified in domhgr. Such modifications break the1030 !! relationship between e12t and e1u*e2u etc.1031 !! Recompute some scale factors ignoring the modified metric.1032 !!----------------------------------------------------------------------1033 !! * Arguments1034 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated1035 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e31036 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors1037 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'1038 !! * Local declarations1039 INTEGER :: ji, jj, jk ! dummy loop indices1040 INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices1041 INTEGER :: isrow ! index for ORCA1 starting row1042 !! acc1043 !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for1044 !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations1045 !!1046 ! ! =====================1047 IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN ! ORCA R2 configuration1048 ! ! =====================1049 !! acc1050 IF( nn_cla == 0 ) THEN1051 !1052 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified)1053 ij0 = 102 ; ij1 = 1021054 DO jk = 1, jpkm11055 DO jj = mj0(ij0), mj1(ij1)1056 DO ji = mi0(ii0), mi1(ii1)1057 SELECT CASE ( pout )1058 CASE( 'U' )1059 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1060 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1061 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1062 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1063 CASE( 'F' )1064 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1065 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1066 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1067 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1068 END SELECT1069 END DO1070 END DO1071 END DO1072 !1073 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified)1074 ij0 = 88 ; ij1 = 881075 DO jk = 1, jpkm11076 DO jj = mj0(ij0), mj1(ij1)1077 DO ji = mi0(ii0), mi1(ii1)1078 SELECT CASE ( pout )1079 CASE( 'U' )1080 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1081 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1082 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1083 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1084 CASE( 'V' )1085 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1086 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1087 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1088 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1089 CASE( 'F' )1090 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1091 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1092 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1093 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1094 END SELECT1095 END DO1096 END DO1097 END DO1098 ENDIF1099 1100 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified)1101 ij0 = 116 ; ij1 = 1161102 DO jk = 1, jpkm11103 DO jj = mj0(ij0), mj1(ij1)1104 DO ji = mi0(ii0), mi1(ii1)1105 SELECT CASE ( pout )1106 CASE( 'U' )1107 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1108 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1109 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1110 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1111 CASE( 'F' )1112 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1113 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1114 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1115 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1116 END SELECT1117 END DO1118 END DO1119 END DO1120 ENDIF1121 !1122 ! ! =====================1123 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration1124 ! ! =====================1125 ! This dirty section will be suppressed by simplification process:1126 ! all this will come back in input files1127 ! Currently these hard-wired indices relate to configuration with1128 ! extend grid (jpjglo=332)1129 ! which had a grid-size of 362x292.1130 isrow = 332 - jpjglo1131 !1132 ii0 = 282 ; ii1 = 283 ! Gibraltar Strait (e2u was modified)1133 ij0 = 241 - isrow ; ij1 = 241 - isrow1134 DO jk = 1, jpkm11135 DO jj = mj0(ij0), mj1(ij1)1136 DO ji = mi0(ii0), mi1(ii1)1137 SELECT CASE ( pout )1138 CASE( 'U' )1139 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1140 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1141 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1142 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1143 CASE( 'F' )1144 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1145 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1146 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1147 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1148 END SELECT1149 END DO1150 END DO1151 END DO1152 !1153 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified)1154 ij0 = 248 - isrow ; ij1 = 248 - isrow1155 DO jk = 1, jpkm11156 DO jj = mj0(ij0), mj1(ij1)1157 DO ji = mi0(ii0), mi1(ii1)1158 SELECT CASE ( pout )1159 CASE( 'U' )1160 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1161 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1162 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1163 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1164 CASE( 'F' )1165 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1166 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1167 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1168 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1169 END SELECT1170 END DO1171 END DO1172 END DO1173 !1174 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified)1175 ij0 = 164 - isrow ; ij1 = 165 - isrow1176 DO jk = 1, jpkm11177 DO jj = mj0(ij0), mj1(ij1)1178 DO ji = mi0(ii0), mi1(ii1)1179 SELECT CASE ( pout )1180 CASE( 'V' )1181 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1182 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1183 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1184 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1185 END SELECT1186 END DO1187 END DO1188 END DO1189 !1190 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on]1191 ij0 = 164 - isrow ; ij1 = 165 - isrow1192 DO jk = 1, jpkm11193 DO jj = mj0(ij0), mj1(ij1)1194 DO ji = mi0(ii0), mi1(ii1)1195 SELECT CASE ( pout )1196 CASE( 'V' )1197 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1198 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1199 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1200 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1201 END SELECT1202 END DO1203 END DO1204 END DO1205 !1206 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified)1207 ij0 = 164 - isrow ; ij1 = 165 - isrow1208 DO jk = 1, jpkm11209 DO jj = mj0(ij0), mj1(ij1)1210 DO ji = mi0(ii0), mi1(ii1)1211 SELECT CASE ( pout )1212 CASE( 'V' )1213 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1214 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1215 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1216 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1217 END SELECT1218 END DO1219 END DO1220 END DO1221 !1222 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified)1223 ij0 = 164 - isrow ; ij1 = 165 - isrow1224 DO jk = 1, jpkm11225 DO jj = mj0(ij0), mj1(ij1)1226 DO ji = mi0(ii0), mi1(ii1)1227 SELECT CASE ( pout )1228 CASE( 'V' )1229 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1230 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1231 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1232 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1233 END SELECT1234 END DO1235 END DO1236 END DO1237 !1238 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified)1239 ij0 = 181 - isrow ; ij1 = 182 - isrow1240 DO jk = 1, jpkm11241 DO jj = mj0(ij0), mj1(ij1)1242 DO ji = mi0(ii0), mi1(ii1)1243 SELECT CASE ( pout )1244 CASE( 'V' )1245 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1246 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1247 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1248 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1249 END SELECT1250 END DO1251 END DO1252 END DO1253 !1254 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified)1255 ij0 = 181 - isrow ; ij1 = 182 - isrow1256 DO jk = 1, jpkm11257 DO jj = mj0(ij0), mj1(ij1)1258 DO ji = mi0(ii0), mi1(ii1)1259 SELECT CASE ( pout )1260 CASE( 'V' )1261 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1262 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1263 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1264 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1265 END SELECT1266 END DO1267 END DO1268 END DO1269 ENDIF1270 ! ! =====================1271 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration1272 ! ! =====================1273 !1274 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified)1275 ij0 = 327 ; ij1 = 3271276 DO jk = 1, jpkm11277 DO jj = mj0(ij0), mj1(ij1)1278 DO ji = mi0(ii0), mi1(ii1)1279 SELECT CASE ( pout )1280 CASE( 'U' )1281 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1282 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1283 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1284 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1285 CASE( 'F' )1286 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1287 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1288 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1289 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1290 END SELECT1291 END DO1292 END DO1293 END DO1294 !1295 ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified)1296 ij0 = 343 ; ij1 = 3431297 DO jk = 1, jpkm11298 DO jj = mj0(ij0), mj1(ij1)1299 DO ji = mi0(ii0), mi1(ii1)1300 SELECT CASE ( pout )1301 CASE( 'U' )1302 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1303 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1304 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1305 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1306 CASE( 'F' )1307 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1308 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1309 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1310 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1311 END SELECT1312 END DO1313 END DO1314 END DO1315 !1316 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified)1317 ij0 = 232 ; ij1 = 2321318 DO jk = 1, jpkm11319 DO jj = mj0(ij0), mj1(ij1)1320 DO ji = mi0(ii0), mi1(ii1)1321 SELECT CASE ( pout )1322 CASE( 'U' )1323 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1324 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1325 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1326 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1327 CASE( 'F' )1328 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1329 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1330 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1331 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1332 END SELECT1333 END DO1334 END DO1335 END DO1336 !1337 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified)1338 ij0 = 232 ; ij1 = 2321339 DO jk = 1, jpkm11340 DO jj = mj0(ij0), mj1(ij1)1341 DO ji = mi0(ii0), mi1(ii1)1342 SELECT CASE ( pout )1343 CASE( 'U' )1344 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1345 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1346 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1347 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1348 CASE( 'F' )1349 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1350 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1351 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1352 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1353 END SELECT1354 END DO1355 END DO1356 END DO1357 !1358 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified)1359 ij0 = 270 ; ij1 = 2701360 DO jk = 1, jpkm11361 DO jj = mj0(ij0), mj1(ij1)1362 DO ji = mi0(ii0), mi1(ii1)1363 SELECT CASE ( pout )1364 CASE( 'U' )1365 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) &1366 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &1367 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &1368 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk)1369 CASE( 'F' )1370 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) &1371 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) &1372 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &1373 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk)1374 END SELECT1375 END DO1376 END DO1377 END DO1378 !1379 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified)1380 ij0 = 232 ; ij1 = 2331381 DO jk = 1, jpkm11382 DO jj = mj0(ij0), mj1(ij1)1383 DO ji = mi0(ii0), mi1(ii1)1384 SELECT CASE ( pout )1385 CASE( 'V' )1386 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1387 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1388 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1389 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1390 END SELECT1391 END DO1392 END DO1393 END DO1394 !1395 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified)1396 ij0 = 276 ; ij1 = 2761397 DO jk = 1, jpkm11398 DO jj = mj0(ij0), mj1(ij1)1399 DO ji = mi0(ii0), mi1(ii1)1400 SELECT CASE ( pout )1401 CASE( 'V' )1402 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) &1403 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &1404 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &1405 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk)1406 END SELECT1407 END DO1408 END DO1409 END DO1410 ENDIF1411 END SUBROUTINE dom_vvl_orca_fix1412 1020 1413 1021 !!====================================================================== 1414 1022 END MODULE domvvl 1415 1416 1417
Note: See TracChangeset
for help on using the changeset viewer.
