- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF
- Files:
-
- 10 edited
-
zdf_oce.F90 (modified) (2 diffs)
-
zdfbfr.F90 (modified) (13 diffs)
-
zdfddm.F90 (modified) (12 diffs)
-
zdfevd.F90 (modified) (2 diffs)
-
zdfgls.F90 (modified) (1 diff)
-
zdfini.F90 (modified) (2 diffs)
-
zdfkpp.F90 (modified) (12 diffs)
-
zdfmxl.F90 (modified) (6 diffs)
-
zdftke.F90 (modified) (15 diffs)
-
zdftmx.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
r4147 r5208 40 40 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 41 41 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: bfrua, bfrva !: Bottom friction coefficients set in zdfbfr 42 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: tfrua, tfrva !: top friction coefficients set in zdfbfr 42 43 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 43 44 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm , avt !: vertical viscosity & diffusivity coef at w-pt [m2/s] … … 57 58 ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) , & 58 59 & avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) , & 60 & tfrua(jpi, jpj), tfrva(jpi, jpj) , & 59 61 & avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk) , & 60 62 & avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk) , STAT = zdf_oce_alloc ) -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r4624 r5208 41 41 REAL(wp), PUBLIC :: rn_bfrien ! local factor to enhance coefficient bfri (PUBLIC for TAM) 42 42 LOGICAL , PUBLIC :: ln_bfr2d ! logical switch for 2D enhancement (PUBLIC for TAM) 43 REAL(wp), PUBLIC :: rn_tfri1 ! top drag coefficient (linear case) (PUBLIC for TAM) 44 REAL(wp), PUBLIC :: rn_tfri2 ! top drag coefficient (non linear case) (PUBLIC for TAM) 45 REAL(wp), PUBLIC :: rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM) 46 REAL(wp), PUBLIC :: rn_tfeb2 ! background top turbulent kinetic energy [m2/s2] (PUBLIC for TAM) 47 REAL(wp), PUBLIC :: rn_tfrien ! local factor to enhance coefficient tfri (PUBLIC for TAM) 48 LOGICAL , PUBLIC :: ln_tfr2d ! logical switch for 2D enhancement (PUBLIC for TAM) 49 43 50 LOGICAL , PUBLIC :: ln_loglayer ! switch for log layer bfr coeff. (PUBLIC for TAM) 44 51 REAL(wp), PUBLIC :: rn_bfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 52 REAL(wp), PUBLIC :: rn_tfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 45 53 LOGICAL , PUBLIC :: ln_bfrimp ! logical switch for implicit bottom friction 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d ! 2D bottomdrag coefficient (PUBLIC for TAM)54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d, tfrcoef2d ! 2D bottom/top drag coefficient (PUBLIC for TAM) 47 55 48 56 !! * Substitutions … … 60 68 !! *** FUNCTION zdf_bfr_alloc *** 61 69 !!---------------------------------------------------------------------- 62 ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )70 ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 63 71 ! 64 72 IF( lk_mpp ) CALL mpp_sum ( zdf_bfr_alloc ) … … 88 96 INTEGER :: ikbt, ikbu, ikbv ! local integers 89 97 REAL(wp) :: zvu, zuv, zecu, zecv, ztmp ! temporary scalars 90 REAL(wp), POINTER, DIMENSION(:,:) :: zbfrt 98 REAL(wp), POINTER, DIMENSION(:,:) :: zbfrt, ztfrt 91 99 !!---------------------------------------------------------------------- 92 100 ! … … 101 109 IF( nn_bfr == 2 ) THEN ! quadratic bottom friction only 102 110 ! 103 CALL wrk_alloc( jpi, jpj, zbfrt )111 CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt ) 104 112 105 113 IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient 106 114 107 # if defined key_vectopt_loop108 DO jj = 1, 1109 !CDIR NOVERRCHK110 DO ji = 1, jpij ! vector opt. (forced unrolling)111 # else112 !CDIR NOVERRCHK113 115 DO jj = 1, jpj 114 !CDIR NOVERRCHK115 116 DO ji = 1, jpi 116 # endif117 117 ikbt = mbkt(ji,jj) 118 ! JC: possible WAD implementation should modify line below if layers vanish118 !! JC: possible WAD implementation should modify line below if layers vanish 119 119 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 120 120 zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 121 121 zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 122 ! (ISF) 123 ikbt = mikt(ji,jj) 124 ! JC: possible WAD implementation should modify line below if layers vanish 125 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 126 ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 127 ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 128 122 129 END DO 123 130 END DO … … 125 132 ELSE 126 133 zbfrt(:,:) = bfrcoef2d(:,:) 127 ENDIF 128 129 # if defined key_vectopt_loop 130 DO jj = 1, 1 131 !CDIR NOVERRCHK 132 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 133 # else 134 !CDIR NOVERRCHK 134 ztfrt(:,:) = tfrcoef2d(:,:) 135 ENDIF 136 135 137 DO jj = 2, jpjm1 136 !CDIR NOVERRCHK137 138 DO ji = 2, jpim1 138 # endif139 139 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 140 140 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) … … 150 150 bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) * zecu 151 151 bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) * zecv 152 ! 153 ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 154 IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 155 bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) & 156 & + ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) ) & 157 & * zecu * (1._wp - umask(ji,jj,1)) 158 END IF 159 IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 160 bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) & 161 & + ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) ) & 162 & * zecv * (1._wp - vmask(ji,jj,1)) 163 END IF 164 ! (ISF) ======================================================================== 165 ikbu = miku(ji,jj) ! ocean bottom level at u- and v-points 166 ikbv = mikv(ji,jj) ! (deepest ocean u- and v-points) 167 ! 168 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & 169 & + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) ) 170 zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) & 171 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 172 ! 173 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 174 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 175 ! 176 tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) * zecu * (1._wp - umask(ji,jj,1)) 177 tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 178 ! (ISF) END ==================================================================== 179 ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 180 IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 181 tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) & 182 & + ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) ) & 183 & * zecu * (1._wp - umask(ji,jj,1)) 184 END IF 185 IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 186 tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) & 187 & + ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) ) & 188 & * zecv * (1._wp - vmask(ji,jj,1)) 189 END IF 152 190 END DO 153 191 END DO 154 155 192 ! 156 193 CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition … … 158 195 IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, & 159 196 & tab2d_2=bfrva, clinfo2= ' v: ', mask2=vmask,ovlap=1 ) 160 CALL wrk_dealloc( jpi,jpj,zbfrt )197 CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt ) 161 198 ENDIF 162 199 ! … … 183 220 INTEGER :: ios 184 221 REAL(wp) :: zminbfr, zmaxbfr ! temporary scalars 222 REAL(wp) :: zmintfr, zmaxtfr ! temporary scalars 185 223 REAL(wp) :: ztmp, zfru, zfrv ! - - 186 224 !! 187 225 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 188 & rn_bfrien, ln_bfrimp, ln_loglayer 226 & rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, & 227 & rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer 189 228 !!---------------------------------------------------------------------- 190 229 ! … … 215 254 bfrua(:,:) = 0.e0 216 255 bfrva(:,:) = 0.e0 256 tfrua(:,:) = 0.e0 257 tfrva(:,:) = 0.e0 217 258 ! 218 259 CASE( 1 ) 219 260 IF(lwp) WRITE(numout,*) ' linear botton friction' 220 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri1 = ', rn_bfri1261 IF(lwp) WRITE(numout,*) ' bottom friction coef. rn_bfri1 = ', rn_bfri1 221 262 IF( ln_bfr2d ) THEN 222 263 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d 223 264 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 265 ENDIF 266 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_bfri1 267 IF( ln_tfr2d ) THEN 268 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 269 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 224 270 ENDIF 225 271 ! … … 252 298 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 253 299 ENDIF 300 IF(lwp) WRITE(numout,*) ' quadratic top friction' 301 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_tfri2 302 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 303 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 304 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 305 IF(lwp) WRITE(numout,*) ' bottom roughness rn_tfrz0 [m] = ', rn_tfrz0 306 IF( rn_tfrz0<=0.e0 ) THEN 307 WRITE(ctmp1,*) ' bottom roughness must be strictly positive' 308 CALL ctl_stop( ctmp1 ) 309 ENDIF 310 IF( ln_tfr2d ) THEN 311 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 312 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 313 ENDIF 254 314 ! 255 315 IF(ln_bfr2d) THEN … … 265 325 ! 266 326 IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 267 # if defined key_vectopt_loop268 DO jj = 1, 1269 !CDIR NOVERRCHK270 DO ji = 1, jpij ! vector opt. (forced unrolling)271 # else272 !CDIR NOVERRCHK273 327 DO jj = 1, jpj 274 !CDIR NOVERRCHK275 328 DO ji = 1, jpi 276 # endif277 329 ikbt = mbkt(ji,jj) 278 330 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp … … 308 360 zminbfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient 309 361 zmaxbfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient 310 ! 311 # if defined key_vectopt_loop 312 DO jj = 1, 1 313 !CDIR NOVERRCHK 314 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 315 # else 316 !CDIR NOVERRCHK 362 zmintfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient 363 zmaxtfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient 364 ! 317 365 DO jj = 2, jpjm1 318 !CDIR NOVERRCHK319 366 DO ji = 2, jpim1 320 # endif321 367 ikbu = mbku(ji,jj) ! deepest ocean level at u- and v-points 322 368 ikbv = mbkv(ji,jj) … … 352 398 WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 353 399 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 400 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 354 401 WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 355 402 ENDIF -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r4624 r5208 6 6 !! History : OPA ! 2000-08 (G. Madec) double diffusive mixing 7 7 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 9 !! 3.6 ! 2013-04 (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_zdfddm || defined key_esopa … … 18 19 USE dom_oce ! ocean space and time domain variables 19 20 USE zdf_oce ! ocean vertical physics variables 21 USE eosbn2 ! equation of state 22 ! 20 23 USE in_out_manager ! I/O manager 21 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 34 37 LOGICAL , PUBLIC, PARAMETER :: lk_zdfddm = .TRUE. !: double diffusive mixing flag 35 38 36 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avs !: salinity vertical diffusivity coeff. at w-point 37 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: rrau !: heat/salt buoyancy flux ratio 38 39 ! !!* Namelist namzdf_ddm : double diffusive mixing * 40 REAL(wp) :: rn_avts ! maximum value of avs for salt fingering 41 REAL(wp) :: rn_hsbfr ! heat/salt buoyancy flux ratio 39 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avs !: salinity vertical diffusivity coeff. at w-point 40 41 ! !!* Namelist namzdf_ddm : double diffusive mixing * 42 REAL(wp) :: rn_avts ! maximum value of avs for salt fingering 43 REAL(wp) :: rn_hsbfr ! heat/salt buoyancy flux ratio 42 44 43 45 !! * Substitutions 46 # include "domzgr_substitute.h90" 44 47 # include "vectopt_loop_substitute.h90" 45 48 !!---------------------------------------------------------------------- 46 !! NEMO/OPA 4.0 , NEMO Consortium (2011)49 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 47 50 !! $Id$ 48 51 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 54 57 !! *** ROUTINE zdf_ddm_alloc *** 55 58 !!---------------------------------------------------------------------- 56 ALLOCATE( avs(jpi,jpj,jpk), rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc ) 57 ! 59 ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc ) 58 60 IF( lk_mpp ) CALL mpp_sum ( zdf_ddm_alloc ) 59 61 IF( zdf_ddm_alloc /= 0 ) CALL ctl_warn('zdf_ddm_alloc: failed to allocate arrays') … … 71 73 !! diffusive mixing (i.e. salt fingering and diffusive layering) 72 74 !! following Merryfield et al. (1999). The rate of double diffusive 73 !! mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S] 74 !! which is computed in rn2.F 75 !! mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]): 75 76 !! * salt fingering (Schmitt 1981): 76 !! for R rau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )77 !! for R rau> 1 and rn2 > 0 : zavfs = O78 !! otherwise : zavft = 0.7 zavs / R rau77 !! for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 ) 78 !! for R > 1 and rn2 > 0 : zavfs = O 79 !! otherwise : zavft = 0.7 zavs / R 79 80 !! * diffusive layering (Federov 1988): 80 !! for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 81 !! * exp( 4.6 exp(-0.54 (1/Rrau-1) ) ) 81 !! for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) ) 82 82 !! otherwise : zavdt = 0 83 !! for .5 < R rau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau-0.85)84 !! for 0 < R rau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau83 !! for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85) 84 !! for 0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R 85 85 !! otherwise : zavds = 0 86 86 !! * update the eddy diffusivity: … … 96 96 ! 97 97 INTEGER :: ji, jj , jk ! dummy loop indices 98 REAL(wp) :: zinr, zrr ! temporary scalars 99 REAL(wp) :: zavft, zavfs ! - - 100 REAL(wp) :: zavdt, zavds ! - - 101 REAL(wp), POINTER, DIMENSION(:,:) :: zmsks, zmskf, zmskd1, zmskd2, zmskd3 98 REAL(wp) :: zaw, zbw, zrw ! local scalars 99 REAL(wp) :: zdt, zds 100 REAL(wp) :: zinr, zrr ! - - 101 REAL(wp) :: zavft, zavfs ! - - 102 REAL(wp) :: zavdt, zavds ! - - 103 REAL(wp), POINTER, DIMENSION(:,:) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 102 104 !!---------------------------------------------------------------------- 103 105 ! 104 106 IF( nn_timing == 1 ) CALL timing_start('zdf_ddm') 105 107 ! 106 CALL wrk_alloc( jpi,jpj, z msks, zmskf, zmskd1, zmskd2, zmskd3 )107 108 CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 109 ! 108 110 ! ! =============== 109 111 DO jk = 2, jpkm1 ! Horizontal slab … … 111 113 ! Define the mask 112 114 ! --------------- 113 rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) ) ! only retains positive value of rrau 115 DO jj = 1, jpj ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 116 DO ji = 1, jpi 117 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 118 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 119 ! 120 zaw = ( rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw ) & 121 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 122 zbw = ( rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw ) & 123 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 124 ! 125 zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 126 zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 127 IF( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 128 zrau(ji,jj) = MAX( 1.e-20, zdt / zds ) ! only retains positive value of zrau 129 END DO 130 END DO 114 131 115 132 DO jj = 1, jpj ! indicators: … … 119 136 ELSE ; zmsks(ji,jj) = 1._wp 120 137 ENDIF 121 ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere122 IF( rrau(ji,jj,jk) <= 1.) THEN ; zmskf(ji,jj) = 0._wp138 ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere 139 IF( zrau(ji,jj) <= 1. ) THEN ; zmskf(ji,jj) = 0._wp 123 140 ELSE ; zmskf(ji,jj) = 1._wp 124 141 ENDIF 125 142 ! diffusive layering indicators: 126 ! ! mskdl1=1 if 0< rrau<1; 0 elsewhere127 IF( rrau(ji,jj,jk) >= 1.) THEN ; zmskd1(ji,jj) = 0._wp143 ! ! mskdl1=1 if 0< R <1; 0 elsewhere 144 IF( zrau(ji,jj) >= 1. ) THEN ; zmskd1(ji,jj) = 0._wp 128 145 ELSE ; zmskd1(ji,jj) = 1._wp 129 146 ENDIF 130 ! ! mskdl2=1 if 0< rrau<0.5; 0 elsewhere131 IF( rrau(ji,jj,jk) >= 0.5) THEN ; zmskd2(ji,jj) = 0._wp147 ! ! mskdl2=1 if 0< R <0.5; 0 elsewhere 148 IF( zrau(ji,jj) >= 0.5 ) THEN ; zmskd2(ji,jj) = 0._wp 132 149 ELSE ; zmskd2(ji,jj) = 1._wp 133 150 ENDIF 134 ! mskdl3=1 if 0.5< rrau<1; 0 elsewhere135 IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN ; zmskd3(ji,jj) = 0._wp136 ELSE 151 ! mskdl3=1 if 0.5< R <1; 0 elsewhere 152 IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN ; zmskd3(ji,jj) = 0._wp 153 ELSE ; zmskd3(ji,jj) = 1._wp 137 154 ENDIF 138 155 END DO … … 149 166 !CDIR NOVERRCHK 150 167 DO ji = 1, jpi 151 zinr = 1. /rrau(ji,jj,jk)168 zinr = 1._wp / zrau(ji,jj) 152 169 ! salt fingering 153 zrr = rrau(ji,jj,jk)/rn_hsbfr170 zrr = zrau(ji,jj) / rn_hsbfr 154 171 zrr = zrr * zrr 155 172 zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj) … … 157 174 ! diffusive layering 158 175 zavdt = 1.3635e-6 * EXP( 4.6 * EXP( -0.54*(zinr-1.) ) ) * zmsks(ji,jj) * zmskd1(ji,jj) 159 zavds = zavdt * zmsks(ji,jj) * ( ( 1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj) &160 & + 0.15 * rrau(ji,jj,jk) * zmskd2(ji,jj) )176 zavds = zavdt * zmsks(ji,jj) * ( ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj) & 177 & + 0.15 * zrau(ji,jj) * zmskd2(ji,jj) ) 161 178 ! add to the eddy viscosity coef. previously computed 162 179 avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds … … 196 213 ENDIF 197 214 ! 198 CALL wrk_dealloc( jpi,jpj, z msks, zmskf, zmskd1, zmskd2, zmskd3 )215 CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 199 216 ! 200 217 IF( nn_timing == 1 ) CALL timing_stop('zdf_ddm') … … 212 229 !! called by zdf_ddm at the first timestep (nit000) 213 230 !!---------------------------------------------------------------------- 231 INTEGER :: ios ! local integer 232 !! 214 233 NAMELIST/namzdf_ddm/ rn_avts, rn_hsbfr 215 INTEGER :: ios ! Local integer output status for namelist read216 234 !!---------------------------------------------------------------------- 217 235 ! -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r3294 r5208 78 78 ! 79 79 DO jk = 1, jpkm1 80 #if defined key_vectopt_loop81 DO jj = 1, 1 ! big loop forced82 DO ji = jpi+2, jpij83 #else84 80 DO jj = 2, jpj ! no vector opt. 85 81 DO ji = 2, jpi 86 #endif87 82 #if defined key_zdfkpp 88 83 ! no evd mixing in the boundary layer with KPP … … 110 105 DO jk = 1, jpkm1 111 106 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 112 #if defined key_vectopt_loop113 DO jj = 1, 1 ! big loop forced114 DO ji = 1, jpij115 #else116 107 DO jj = 1, jpj ! loop over the whole domain (no lbc_lnk call) 117 108 DO ji = 1, jpi 118 #endif119 109 #if defined key_zdfkpp 120 110 ! no evd mixing in the boundary layer with KPP -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r4624 r5208 1258 1258 en (:,:,:) = rn_emin 1259 1259 mxln(:,:,:) = 0.001 1260 avt_k (:,:,:) = avt (:,:,:) 1261 avm_k (:,:,:) = avm (:,:,:) 1262 avmu_k(:,:,:) = avmu(:,:,:) 1263 avmv_k(:,:,:) = avmv(:,:,:) 1260 1264 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_gls( jit ) ; END DO 1261 1265 ENDIF -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
r4677 r5208 14 14 !!---------------------------------------------------------------------- 15 15 USE par_oce ! mesh and scale factors 16 USE sbc_oce ! surface module (only for nn_isf in the option compatibility test) 16 17 USE ldftra_oce ! ocean active tracers: lateral physics 17 18 USE ldfdyn_oce ! ocean dynamics lateral physics … … 117 118 IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa ) & 118 119 & CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 120 IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 ) & 121 & CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 119 122 ! 120 123 ! ! ... Convection -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
r4624 r5208 26 26 USE phycst ! physical constants 27 27 USE eosbn2 ! equation of state 28 USE zdfddm ! double diffusion mixing 28 USE zdfddm ! double diffusion mixing (avs array) 29 USE lib_mpp ! MPP library 30 USE trd_oce ! ocean trends definition 31 USE trdtra ! tracers trends 32 ! 29 33 USE in_out_manager ! I/O manager 30 USE lib_mpp ! MPP library31 USE wrk_nemo ! work arrays32 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 35 USE prtctl ! Print control 34 USE trdmod_oce ! ocean trends definition 35 USE trdtra ! tracers trends 36 USE wrk_nemo ! work arrays 36 37 USE timing ! Timing 37 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 38 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 38 39 39 40 IMPLICIT NONE … … 246 247 REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t 247 248 #if defined key_zdfddm 248 REAL(wp) :: zrrau, zds, zavdds, zavddt,zinr ! double diffusion mixing 249 REAL(wp), POINTER, DIMENSION(:,:) :: zdifs 249 REAL(wp) :: zrw, zkm1s ! local scalars 250 REAL(wp) :: zrrau, zdt, zds, zavdds, zavddt, zinr ! double diffusion mixing 251 REAL(wp), POINTER, DIMENSION(:,:) :: zdifs 250 252 REAL(wp), POINTER, DIMENSION(:) :: za2s, za3s, zkmps 251 REAL(wp) :: zkm1s252 253 REAL(wp), POINTER, DIMENSION(:,:) :: zblcs 253 254 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdiffus … … 274 275 #endif 275 276 276 zviscos(:,:,:) = 0. 277 zblcm (:,: ) = 0. 278 zdiffut(:,:,:) = 0. 279 zblct (:,: ) = 0. 277 zviscos(:,:,:) = 0._wp 278 zblcm (:,: ) = 0._wp 279 zdiffut(:,:,:) = 0._wp 280 zblct (:,: ) = 0._wp 280 281 #if defined key_zdfddm 281 zdiffus(:,:,:) = 0. 282 zblcs (:,: ) = 0. 283 #endif 284 ghats(:,:,:) = 0. 285 286 zBo (:,:) = 0. 287 zBosol(:,:) = 0. 288 zustar(:,:) = 0. 289 290 282 zdiffus(:,:,:) = 0._wp 283 zblcs (:,: ) = 0._wp 284 #endif 285 ghats (:,:,:) = 0._wp 286 zBo (:,: ) = 0._wp 287 zBosol (:,: ) = 0._wp 288 zustar (:,: ) = 0._wp 289 ! 291 290 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 292 291 ! I. Interior diffusivity and viscosity at w points ( T interfaces) … … 332 331 avt (ji,jj,jk) = avt (ji,jj,jk) + rn_difri * zfri 333 332 ENDIF 333 ! 334 334 #if defined key_zdfddm 335 avs (ji,jj,jk) = avt (ji,jj,jk)335 ! 336 336 ! Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90 337 ! ------------------------------------------------------------------ 338 ! only retains positive value of rrau 339 zrrau = MAX( rrau(ji,jj,jk), epsln ) 340 zds = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) 341 IF( zrrau > 1. .AND. zds > 0.) THEN 342 ! 343 ! Salt fingering case. 344 !--------------------- 345 ! Compute interior diffusivity for double diffusive mixing of 346 ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 347 ! After that set interior diffusivity for double diffusive mixing 348 ! of temperature 337 ! ------------------------- 338 avs (ji,jj,jk) = avt (ji,jj,jk) 339 340 ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 341 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 342 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 343 ! 344 zaw = ( rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw ) * tmask(ji,jj,jk) 345 zbw = ( rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw ) * tmask(ji,jj,jk) 346 ! 347 zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 348 zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 349 IF( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 350 zrrau = MAX( epsln , zdt / zds ) ! only retains positive value of zrau 351 ! 352 IF( zrrau > 1. .AND. zds > 0.) THEN ! Salt fingering case. 353 ! !--------------------- 354 ! Compute interior diffusivity for double diffusive mixing of salinity. 355 ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 356 ! After that set interior diffusivity for double diffusive mixing of temperature 349 357 zavdds = MIN( zrrau, Rrho0 ) 350 358 zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 ) … … 353 361 zavdds = difssf * zavdds 354 362 zavddt = 0.7 * zavdds 355 ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN356 363 ! 357 ! Diffusive convection case. 358 !--------------------------- 359 ! Compute interior diffusivity for double diffusive mixing of 360 ! temperature (Marmorino and Caldwell, 1976); 364 ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN ! Diffusive convection case. 365 ! !--------------------------- 366 ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976); 361 367 ! Compute interior diffusivity for double diffusive mixing of salinity 362 368 zinr = 1. / zrrau 363 369 zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 364 370 zavddt = difsdc * zavddt 365 IF( zrrau < 0.5) THEN 366 zavdds = zavddt * 0.15 * zrrau 367 ELSE 368 zavdds = zavddt * (1.85 * zrrau - 0.85 ) 371 IF( zrrau < 0.5) THEN ; zavdds = zavddt * 0.15 * zrrau 372 ELSE ; zavdds = zavddt * (1.85 * zrrau - 0.85 ) 369 373 ENDIF 370 374 ELSE … … 385 389 !--------------------------------------------------------------------- 386 390 DO jj = 2, jpjm1 387 DO ji = fs_2, fs_jpim1 388 IF( nn_eos < 1) THEN 389 zt = tsn(ji,jj,1,jp_tem) 390 zs = tsn(ji,jj,1,jp_sal) - 35.0 391 zh = fsdept(ji,jj,1) 392 ! potential volumic mass 393 zrhos = rhop(ji,jj,1) 394 zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta 395 & - 0.203814e-03 ) * zt & 396 & + 0.170907e-01 ) * zt & 397 & + 0.665157e-01 & 398 & + ( - 0.678662e-05 * zs & 399 & - 0.846960e-04 * zt + 0.378110e-02 ) * zs & 400 & + ( ( - 0.302285e-13 * zh & 401 & - 0.251520e-11 * zs & 402 & + 0.512857e-12 * zt * zt ) * zh & 403 & - 0.164759e-06 * zs & 404 & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & 405 & + 0.380374e-04 ) * zh 406 407 zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta 408 & - 0.301985e-05 ) * zt & 409 & + 0.785567e-03 & 410 & + ( 0.515032e-08 * zs & 411 & + 0.788212e-08 * zt - 0.356603e-06 ) * zs & 412 & +( ( 0.121551e-17 * zh & 413 & - 0.602281e-15 * zs & 414 & - 0.175379e-14 * zt + 0.176621e-12 ) * zh & 415 & + 0.408195e-10 * zs & 416 & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & 417 & - 0.121555e-07 ) * zh 418 419 zthermal = zbeta * zalbet / ( rcp * zrhos + epsln ) 420 zhalin = zbeta * tsn(ji,jj,1,jp_sal) * rcs 421 ELSE 422 zrhos = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) ) 423 zthermal = rn_alpha / ( rcp * zrhos + epsln ) 424 zhalin = rn_beta * tsn(ji,jj,1,jp_sal) * rcs 425 zbeta = rn_beta 426 ENDIF 391 DO ji = fs_2, fs_jpim1 392 zrhos = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1) 393 zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln ) 394 zbeta = rab_n(ji,jj,1,jp_sal) 395 zhalin = zbeta * tsn(ji,jj,1,jp_sal) * rcs 396 ! 427 397 ! Radiative surface buoyancy force 428 398 zBosol(ji,jj) = grav * zthermal * qsr(ji,jj) … … 435 405 ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) & 436 406 & + sfx(ji,jj) ) * rcs * tmask(ji,jj,1) 437 END DO438 END DO407 END DO 408 END DO 439 409 440 410 zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) … … 447 417 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 448 418 zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos + epsln ) ) 449 END DO450 END DO419 END DO 420 END DO 451 421 452 422 !CDIR NOVERRCHK … … 1270 1240 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 1271 1241 !!bug gm jpttdzdf ==> jpttkpp 1272 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_zdf, ztrdt )1273 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_zdf, ztrds )1242 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 1243 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 1274 1244 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1275 1245 ENDIF … … 1340 1310 IF( l_trdtrc ) THEN ! save the non-local tracer flux trends for diagnostic 1341 1311 ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:) 1342 CALL trd_tra( kt, 'TRC', jn, jptra_ trd_zdf, ztrtrd(:,:,:) )1312 CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:) ) 1343 1313 ENDIF 1344 1314 ! … … 1375 1345 !!---------------------------------------------------------------------- 1376 1346 INTEGER :: ji, jj, jk ! dummy loop indices 1347 INTEGER :: ios ! local integer 1377 1348 #if ! defined key_kppcustom 1378 1349 INTEGER :: jm ! dummy loop indices … … 1382 1353 REAL(wp) :: zustar, zucube, zustvk, zeta, zehat ! tempory scalars 1383 1354 #endif 1384 INTEGER :: ios ! Local integer output status for namelist read1385 1355 REAL(wp) :: zhbf ! tempory scalars 1386 1356 LOGICAL :: ll_kppcustom ! 1st ocean level taken as surface layer -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r4245 r5208 6 6 !! History : 1.0 ! 2003-08 (G. Madec) original code 7 7 !! 3.2 ! 2009-07 (S. Masson, G. Madec) IOM + merge of DO-loop 8 !! 3.7 ! 2012-03 (G. Madec) make public the density criteria for trdmxl 9 !! - ! 2014-02 (F. Roquet) mixed layer depth calculated using N2 instead of rhop 8 10 !!---------------------------------------------------------------------- 9 11 !! zdf_mxl : Compute the turbocline and mixed layer depths. … … 14 16 USE in_out_manager ! I/O manager 15 17 USE prtctl ! Print control 18 USE phycst ! physical constants 16 19 USE iom ! I/O library 17 20 USE lib_mpp ! MPP library … … 25 28 PUBLIC zdf_mxl ! called by step.F90 26 29 27 REAL(wp), PUBLIC :: rho_c = 0.01_wp ! density criterion for mixed layer depth28 REAL(wp), PUBLIC :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth29 30 30 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by TOP) 31 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld !: mixing layer depth (turbocline) [m] 32 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 33 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m] 34 35 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 36 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 34 37 35 38 !! * Substitutions … … 70 73 !! eddy diffusivity coefficient (resulting from the vertical physics 71 74 !! alone, not the isopycnal part, see trazdf.F) fall below a given 72 !! value defined locally (avt_c here taken equal to 5 cm/s2 )75 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default) 73 76 !! 74 77 !! ** Action : nmln, hmld, hmlp, hmlpt 75 78 !!---------------------------------------------------------------------- 76 79 INTEGER, INTENT(in) :: kt ! ocean time-step index 77 !! 78 INTEGER :: ji, jj, jk ! dummy loop indices 79 INTEGER :: iikn, iiki ! temporary integer within a do loop 80 INTEGER, POINTER, DIMENSION(:,:) :: imld ! temporary workspace 80 ! 81 INTEGER :: ji, jj, jk ! dummy loop indices 82 INTEGER :: iikn, iiki, ikt, imkt ! local integer 83 REAL(wp) :: zN2_c ! local scalar 84 INTEGER, POINTER, DIMENSION(:,:) :: imld ! 2D workspace 81 85 !!---------------------------------------------------------------------- 82 86 ! … … 94 98 95 99 ! w-level of the mixing and mixed layers 96 nmln(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 97 imld(:,:) = mbkt(:,:) + 1 98 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 100 nmln(:,:) = nlb10 ! Initialization to the number of w ocean point 101 hmlp(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2 102 zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria 103 DO jk = nlb10, jpkm1 104 DO jj = 1, jpj ! Mixed layer level: w-level 105 DO ji = 1, jpi 106 ikt = mbkt(ji,jj) 107 hmlp(ji,jj) = hmlp(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * fse3w(ji,jj,jk) 108 IF( hmlp(ji,jj) < zN2_c ) nmln(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 109 END DO 110 END DO 111 END DO 112 ! 113 ! w-level of the turbocline 114 imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 115 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 99 116 DO jj = 1, jpj 100 117 DO ji = 1, jpi 101 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c ) nmln(ji,jj) = jk ! Mixed layer102 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = jk! Turbocline118 imkt = mikt(ji,jj) 119 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = MAX( imkt, jk ) ! Turbocline 103 120 END DO 104 121 END DO … … 109 126 iiki = imld(ji,jj) 110 127 iikn = nmln(ji,jj) 111 hmld (ji,jj) = fsdepw(ji,jj,iiki ) * tmask(ji,jj,1) ! Turbocline depth 112 hmlp (ji,jj) = fsdepw(ji,jj,iikn ) * tmask(ji,jj,1) ! Mixed layer depth 113 hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) ! depth of the last T-point inside the mixed layer 128 imkt = mikt(ji,jj) 129 hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Turbocline depth 130 hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj) ! Mixed layer depth 131 hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 114 132 END DO 115 133 END DO -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r5208 241 241 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 243 IF (mikt(ji,jj) .GT. 1) THEN 244 en(ji,jj,mikt(ji,jj))=rn_emin 245 ELSE 246 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 247 END IF 244 248 END DO 245 249 END DO … … 291 295 END DO 292 296 ! ! finite LC depth 293 # if defined key_vectopt_loop294 DO jj = 1, 1295 DO ji = 1, jpij ! vector opt. (forced unrolling)296 # else297 297 DO jj = 1, jpj 298 298 DO ji = 1, jpi 299 # endif300 299 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 301 300 END DO 302 301 END DO 303 302 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 304 !CDIR NOVERRCHK305 303 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 306 !CDIR NOVERRCHK307 304 DO jj = 2, jpjm1 308 !CDIR NOVERRCHK309 305 DO ji = fs_2, fs_jpim1 ! vector opt. 310 306 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 342 338 END DO 343 339 ! 344 DO j k = 2, jpkm1 !* Matrix and right hand side in en345 DO j j = 2, jpjm1346 DO j i = fs_2, fs_jpim1 ! vector opt.340 DO jj = 2, jpjm1 341 DO ji = fs_2, fs_jpim1 ! vector opt. 342 DO jk = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en 347 343 zcof = zfact1 * tmask(ji,jj,jk) 348 344 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 360 356 ! ! right hand side in en 361 357 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 362 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) 363 END DO 364 END DO 365 END DO 366 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 367 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 368 DO jj = 2, jpjm1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 358 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 359 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 360 END DO 361 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 362 DO jk = mikt(ji,jj)+2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 370 363 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 371 364 END DO 372 END DO 373 END DO 374 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 375 DO ji = fs_2, fs_jpim1 ! vector opt. 376 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 377 END DO 378 END DO 379 DO jk = 3, jpkm1 380 DO jj = 2, jpjm1 381 DO ji = fs_2, fs_jpim1 ! vector opt. 365 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 366 zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj)) ! Surface boudary conditions on tke 367 ! 368 DO jk = mikt(ji,jj)+2, jpkm1 382 369 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 383 370 END DO 384 END DO 385 END DO 386 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 387 DO ji = fs_2, fs_jpim1 ! vector opt. 371 ! 372 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 388 373 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 389 END DO 390 END DO 391 DO jk = jpk-2, 2, -1 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 ! vector opt. 374 ! 375 DO jk = jpk-2, mikt(ji,jj)+1, -1 394 376 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 395 377 END DO 396 END DO 397 END DO 398 DO jk = 2, jpkm1 ! set the minimum value of tke 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 378 ! 379 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 401 380 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 402 381 END DO … … 412 391 DO ji = fs_2, fs_jpim1 ! vector opt. 413 392 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 414 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)393 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 415 394 END DO 416 395 END DO … … 421 400 jk = nmln(ji,jj) 422 401 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 423 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)402 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 424 403 END DO 425 404 END DO … … 433 412 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 434 413 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 435 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress414 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 436 415 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 437 416 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 438 417 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 439 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)418 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 440 419 END DO 441 420 END DO … … 506 485 ! 507 486 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 508 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 509 zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) ) 510 ELSE ! surface set to the minimum value 511 zmxlm(:,:,1) = rn_mxl0 487 DO jj = 2, jpjm1 488 DO ji = fs_2, fs_jpim1 489 IF (mikt(ji,jj) .GT. 1) THEN 490 zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 491 ELSE 492 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 493 zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 494 END IF 495 END DO 496 END DO 497 ELSE 498 DO jj = 2, jpjm1 499 DO ji = fs_2, fs_jpim1 ! surface set to the minimum value 500 zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 501 END DO 502 END DO 512 503 ENDIF 513 504 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 514 505 ! 515 506 !CDIR NOVERRCHK 516 DO j k = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)517 !CDIR NOVERRCHK 518 DO j j = 2, jpjm1519 !CDIR NOVERRCHK520 DO j i = fs_2, fs_jpim1 ! vector opt.507 DO jj = 2, jpjm1 508 !CDIR NOVERRCHK 509 DO ji = fs_2, fs_jpim1 ! vector opt. 510 !CDIR NOVERRCHK 511 DO jk = mikt(ji,jj)+1, jpkm1 ! interior value : l=sqrt(2*e/n^2) 521 512 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 522 513 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 523 514 END DO 515 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 524 516 END DO 525 517 END DO … … 533 525 ! 534 526 CASE ( 0 ) ! bounded by the distance to surface and bottom 535 DO j k = 2, jpkm1536 DO j j = 2, jpjm1537 DO j i = fs_2, fs_jpim1 ! vector opt.538 zemxl = MIN( fsdepw(ji,jj,jk) , zmxlm(ji,jj,jk), &527 DO jj = 2, jpjm1 528 DO ji = fs_2, fs_jpim1 ! vector opt. 529 DO jk = mikt(ji,jj)+1, jpkm1 530 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 539 531 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 540 532 zmxlm(ji,jj,jk) = zemxl … … 545 537 ! 546 538 CASE ( 1 ) ! bounded by the vertical scale factor 547 DO j k = 2, jpkm1548 DO j j = 2, jpjm1549 DO j i = fs_2, fs_jpim1 ! vector opt.539 DO jj = 2, jpjm1 540 DO ji = fs_2, fs_jpim1 ! vector opt. 541 DO jk = mikt(ji,jj)+1, jpkm1 550 542 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 551 543 zmxlm(ji,jj,jk) = zemxl … … 556 548 ! 557 549 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 558 DO j k = 2, jpkm1 ! from the surface to the bottom :559 DO j j = 2, jpjm1560 DO j i = fs_2, fs_jpim1 ! vector opt.550 DO jj = 2, jpjm1 551 DO ji = fs_2, fs_jpim1 ! vector opt. 552 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : 561 553 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 562 554 END DO 563 END DO 564 END DO 565 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 566 DO jj = 2, jpjm1 567 DO ji = fs_2, fs_jpim1 ! vector opt. 555 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 568 556 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 569 557 zmxlm(ji,jj,jk) = zemxl … … 574 562 ! 575 563 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 576 DO j k = 2, jpkm1 ! from the surface to the bottom : lup577 DO j j = 2, jpjm1578 DO j i = fs_2, fs_jpim1 ! vector opt.564 DO jj = 2, jpjm1 565 DO ji = fs_2, fs_jpim1 ! vector opt. 566 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup 579 567 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 580 568 END DO 581 END DO 582 END DO 583 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 569 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 586 570 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 587 571 END DO … … 628 612 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 629 613 ! 630 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 614 DO jj = 2, jpjm1 615 DO ji = fs_2, fs_jpim1 ! vector opt. 616 DO jk = miku(ji,jj)+1, jpkm1 !* vertical eddy viscosity at u- and v-points 617 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 618 END DO 619 DO jk = mikv(ji,jj)+1, jpkm1 620 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 621 END DO 622 END DO 623 END DO 624 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 625 ! 626 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 631 627 DO jj = 2, jpjm1 632 628 DO ji = fs_2, fs_jpim1 ! vector opt. 633 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 634 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 635 END DO 636 END DO 637 END DO 638 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 639 ! 640 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 641 DO jk = 2, jpkm1 642 DO jj = 2, jpjm1 643 DO ji = fs_2, fs_jpim1 ! vector opt. 629 DO jk = mikt(ji,jj)+1, jpkm1 644 630 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 645 631 ! ! shear … … 655 641 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 656 642 # if defined key_c1d 657 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) 658 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)! c1d config. : save Ri643 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 644 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 659 645 # endif 660 646 END DO … … 740 726 ! 741 727 ! !* Check of some namelist values 742 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 743 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 744 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 745 #if ! key_coupled 746 IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 747 #endif 728 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 729 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 730 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 731 IF( nn_etau == 3 .AND. .NOT. lk_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 748 732 749 733 IF( ln_mxl0 ) THEN -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r4624 r5208 126 126 zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column 127 127 DO jk = 2, jpkm1 128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) 128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 129 129 END DO 130 130 … … 135 135 END DO 136 136 137 DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 138 zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s 137 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 138 DO ji = 1, jpi 139 DO jk = mikt(ji,jj)+1, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s 141 END DO 142 END DO 139 143 END DO 140 144 … … 162 166 ! ! Update mixing coefs ! 163 167 ! ! ----------------------- ! 164 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 165 avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk) 166 avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) 167 DO jj = 2, jpjm1 168 DO ji = fs_2, fs_jpim1 ! vector opt. 168 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 169 DO ji = 1, jpi 170 DO jk = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 171 avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 172 avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 173 END DO 174 END DO 175 END DO 176 177 DO jj = 2, jpjm1 178 DO ji = fs_2, fs_jpim1 ! vector opt. 179 DO jk = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 169 180 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * umask(ji,jj,jk) 170 181 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * vmask(ji,jj,jk) … … 237 248 zsum2(:,:) = 0.e0 238 249 DO jk= 2, jpk 239 zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 240 zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) 250 zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 251 zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 241 252 END DO 242 253 DO jj = 1, jpj … … 274 285 zkz(:,:) = 0.e0 ! Associated potential energy consummed over the whole water column 275 286 DO jk = 2, jpkm1 276 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) 287 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 277 288 END DO 278 289 … … 284 295 285 296 DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s 286 zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) ! kz max = 120 cm2/s297 zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1) ! kz max = 120 cm2/s 287 298 END DO 288 299 … … 414 425 ! only the energy available for mixing is taken into account, 415 426 ! (mixing efficiency tidal dissipation efficiency) 416 en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1) 417 427 en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:) 428 429 !============ 430 !TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep? 418 431 ! Vertical structure (az_tmx) 419 432 DO jj = 1, jpj ! part independent of the level 420 433 DO ji = 1, jpi 421 zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean434 zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 422 435 zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) ) 423 436 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj) … … 427 440 DO jj = 1, jpj 428 441 DO ji = 1, jpi 429 az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk) 430 END DO 431 END DO 432 END DO 442 az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk) 443 END DO 444 END DO 445 END DO 446 !=========== 433 447 434 448 IF( nprint == 1 .AND. lwp ) THEN … … 443 457 ztpc = 0.e0 444 458 zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 445 DO j k= 2, jpkm1446 DO j j = 1, jpj447 DO j i = 1, jpi459 DO jj = 1, jpj 460 DO ji = 1, jpi 461 DO jk= mikt(ji,jj)+1, jpkm1 448 462 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 449 463 END DO … … 459 473 zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 ) 460 474 zkz(:,:) = 0.e0 461 DO j k = 2, jpkm1462 DO jj = 1, jpj463 DO ji = 1, jpi475 DO jj = 1, jpj 476 DO ji = 1, jpi 477 DO jk = mikt(ji,jj)+1, jpkm1 464 478 zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk) 465 END DO466 END DO479 END DO 480 END DO 467 481 END DO 468 482 ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz … … 484 498 WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 485 499 486 DO jk = 2, jpkm1 487 zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s 500 DO jj = 1, jpj 501 DO ji = 1, jpi 502 DO jk = mikt(ji,jj)+1, jpkm1 503 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s 504 END DO 505 END DO 488 506 END DO 489 507 ztpc = 0.e0
Note: See TracChangeset
for help on using the changeset viewer.
