- Timestamp:
- 2017-05-20T13:49:38+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r7991 r8055 16 16 17 17 !!---------------------------------------------------------------------- 18 !! zdf_ric_init : initialization, namelist read, & parameters control 18 19 !! zdf_ric : update momentum and tracer Kz from the Richardson number 19 !! zdf_ric_init : initialization, namelist read, & parameters control20 !! ric_rst : read/write RIC restart in ocean restart file 20 21 !!---------------------------------------------------------------------- 21 22 USE oce ! ocean dynamics and tracers variables 22 23 USE dom_oce ! ocean space and time domain variables 23 24 USE zdf_oce ! vertical physics: variables 24 USE zdfsh2 ! vertical physics: shear production term of TKE25 25 USE phycst ! physical constants 26 26 USE sbc_oce, ONLY : taum … … 28 28 ! 29 29 USE in_out_manager ! I/O manager 30 USE lbclnk ! ocean lateral boundary condition (or mpp link) 31 USE lib_mpp ! MPP library 30 USE iom ! I/O manager library 32 31 USE timing ! Timing 33 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) … … 37 36 PRIVATE 38 37 39 PUBLIC zdf_ric ! called by step.F90 40 PUBLIC zdf_ric_init ! called by opa.F90 38 PUBLIC zdf_ric ! called by zdfphy.F90 39 PUBLIC ric_rst ! called by zdfphy.F90 40 PUBLIC zdf_ric_init ! called by nemogcm.F90 41 41 42 42 ! !!* Namelist namzdf_ric : Richardson number dependent Kz * … … 52 52 53 53 !! * Substitutions 54 # include " vectopt_loop_substitute.h90"54 # include "domain_substitute.h90" 55 55 !!---------------------------------------------------------------------- 56 56 !! NEMO/OPA 4.0 , NEMO Consortium (2011) … … 59 59 !!---------------------------------------------------------------------- 60 60 CONTAINS 61 62 SUBROUTINE zdf_ric( kt )63 !!----------------------------------------------------------------------64 !! *** ROUTINE zdfric ***65 !!66 !! ** Purpose : Compute the before eddy viscosity and diffusivity as67 !! a function of the local richardson number.68 !!69 !! ** Method : Local richardson number dependent formulation of the70 !! vertical eddy viscosity and diffusivity coefficients.71 !! The eddy coefficients are given by:72 !! avm = avm0 + avmb73 !! avt = avm0 / (1 + rn_alp*ri)74 !! with ri = N^2 / dz(u)**275 !! = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]76 !! avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric77 !! Where ri is the before local Richardson number,78 !! rn_avmri is the maximum value reaches by avm and avt79 !! avmb and avtb are the background (or minimum) values80 !! and rn_alp, nn_ric are adjustable parameters.81 !! Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s82 !! avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.83 !! a numerical threshold is impose on the vertical shear (1.e-20)84 !! As second step compute Ekman depth from wind stress forcing85 !! and apply namelist provided vertical coeff within this depth.86 !! The Ekman depth is:87 !! Ustar = SQRT(Taum/rho0)88 !! ekd= rn_ekmfc * Ustar / f089 !! Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation90 !! of the above equation indicates the value is somewhat arbitrary; therefore91 !! we allow the freedom to increase or decrease this value, if the92 !! Ekman depth estimate appears too shallow or too deep, respectively.93 !! Ekd is then limited by rn_mldmin and rn_mldmax provided in the94 !! namelist95 !! N.B. the mask are required for implicit scheme, and surface96 !! and bottom value already set in zdfphy.F9097 !!98 !! ** Action : avm, avt mixing coeff (inner domain values only)99 !!100 !! References : Pacanowski & Philander 1981, JPO, 1441-1451.101 !! PFJ Lermusiaux 2001.102 !!----------------------------------------------------------------------103 INTEGER, INTENT(in) :: kt ! ocean time-step104 !!105 INTEGER :: ji, jj, jk ! dummy loop indices106 LOGICAL :: ll_av_weight = .TRUE. ! local logical107 REAL(wp) :: zcfRi, zav, zustar ! local scalars108 REAL(wp), DIMENSION(jpi,jpj) :: zh_ekm ! 2D workspace109 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsh2 ! 3D -110 !!----------------------------------------------------------------------111 !112 IF( nn_timing == 1 ) CALL timing_start('zdf_ric')113 !114 ! !== avm and avt = F(Richardson number) ==!115 !116 ! !* Shear production at uw- and vw-points (energy conserving form)117 CALL zdf_sh2( zsh2 )118 !119 DO jk = 2, jpkm1 !* Vertical eddy viscosity and diffusivity coefficients120 DO jj = 1, jpjm1121 DO ji = 1, jpim1122 ! ! coefficient = F(richardson number) (avm-weighted Ri)123 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) ) )124 zav = rn_avmri * zcfRi**nn_ric125 !126 ! ! avm and avt coefficients127 avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk)128 avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk)129 END DO130 END DO131 END DO132 !133 !!gm BUG <<<<==== This param can't work at low latitude134 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )135 !136 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==!137 !138 DO jj = 2, jpjm1 !* Ekman depth139 DO ji = 2, jpim1140 zustar = SQRT( taum(ji,jj) * r1_rau0 )141 zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth142 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax ) ) ! set allowed rang143 END DO144 END DO145 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer146 DO jj = 2, jpjm1147 DO ji = 2, jpim1148 IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN149 avm(ji,jj,jk) = MAX( avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk)150 avt(ji,jj,jk) = MAX( avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk)151 ENDIF152 END DO153 END DO154 END DO155 ENDIF156 !157 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric')158 !159 END SUBROUTINE zdf_ric160 161 61 162 62 SUBROUTINE zdf_ric_init … … 205 105 ENDIF 206 106 ! 207 DO jk = 1, jpk ! Initialization of vertical eddy coef. to the background value 208 avt(:,:,jk) = avtb(jk) * tmask(:,:,jk) 209 avm(:,:,jk) = avmb(jk) * umask(:,:,jk) 107 CALL ric_rst( nit000, 'READ' ) !* read or initialize all required files 108 ! 109 END SUBROUTINE zdf_ric_init 110 111 112 SUBROUTINE zdf_ric( ARG_2D, kt, pdept, p_sh2, p_avm, p_avt ) 113 !!---------------------------------------------------------------------- 114 !! *** ROUTINE zdfric *** 115 !! 116 !! ** Purpose : Compute the before eddy viscosity and diffusivity as 117 !! a function of the local richardson number. 118 !! 119 !! ** Method : Local richardson number dependent formulation of the 120 !! vertical eddy viscosity and diffusivity coefficients. 121 !! The eddy coefficients are given by: 122 !! avm = avm0 + avmb 123 !! avt = avm0 / (1 + rn_alp*ri) 124 !! with Ri = N^2 / dz(u)**2 125 !! = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 126 !! avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 127 !! where ri is the before local Richardson number, 128 !! rn_avmri is the maximum value reaches by avm and avt 129 !! and rn_alp, nn_ric are adjustable parameters. 130 !! Typical values : rn_alp=5. and nn_ric=2. 131 !! 132 !! As second step compute Ekman depth from wind stress forcing 133 !! and apply namelist provided vertical coeff within this depth. 134 !! The Ekman depth is: 135 !! Ustar = SQRT(Taum/rho0) 136 !! ekd= rn_ekmfc * Ustar / f0 137 !! Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation 138 !! of the above equation indicates the value is somewhat arbitrary; therefore 139 !! we allow the freedom to increase or decrease this value, if the 140 !! Ekman depth estimate appears too shallow or too deep, respectively. 141 !! Ekd is then limited by rn_mldmin and rn_mldmax provided in the 142 !! namelist 143 !! N.B. the mask are required for implicit scheme, and surface 144 !! and bottom value already set in zdfphy.F90 145 !! 146 !! ** Action : avm, avt mixing coeff (inner domain values only) 147 !! 148 !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 149 !! PFJ Lermusiaux 2001. 150 !!---------------------------------------------------------------------- 151 INTEGER , INTENT(in ) :: ARG_2D ! inner domain start-end i-indices 152 INTEGER , INTENT(in ) :: kt ! ocean time-step 153 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdept ! depth of t-point [m] 154 REAL(wp), DIMENSION(WRK_3D), INTENT(in ) :: p_sh2 ! shear production term 155 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 156 !! 157 INTEGER :: ji, jj, jk ! dummy loop indices 158 REAL(wp) :: zcfRi, zav, zustar, zed ! local scalars 159 REAL(wp), DIMENSION(WRK_2D) :: zh_ekm ! 2D workspace 160 !!---------------------------------------------------------------------- 161 ! 162 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 163 ! 164 ! !== avm and avt = F(Richardson number) ==! 165 DO jk = 2, jpkm1 166 DO jj = k_Jstr, k_Jend 167 DO ji = k_Jstr, k_Iend ! coefficient = F(richardson number) (avm-weighted Ri) 168 zcfRi = 1._wp / ( 1._wp + rn_alp * MAX( 0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) ) ) 169 zav = rn_avmri * zcfRi**nn_ric 170 ! ! avm and avt coefficients 171 p_avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 172 p_avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 173 END DO 174 END DO 210 175 END DO 211 176 ! 212 END SUBROUTINE zdf_ric_init 177 !!gm BUG <<<<==== This param can't work at low latitude 178 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 179 ! 180 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 181 ! 182 DO jj = k_Jstr, k_Jend !* Ekman depth 183 DO ji = k_Jstr, k_Iend 184 zustar = SQRT( taum(ji,jj) * r1_rau0 ) 185 zed = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 186 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zed , rn_mldmax ) ) ! set allowed range 187 END DO 188 END DO 189 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer 190 DO jj = k_Jstr, k_Jend 191 DO ji = k_Jstr, k_Iend 192 IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 193 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) 194 p_avt(ji,jj,jk) = MAX( p_avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk) 195 ENDIF 196 END DO 197 END DO 198 END DO 199 ENDIF 200 ! 201 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') 202 ! 203 END SUBROUTINE zdf_ric 204 205 206 SUBROUTINE ric_rst( kt, cdrw ) 207 !!--------------------------------------------------------------------- 208 !! *** ROUTINE ric_rst *** 209 !! 210 !! ** Purpose : Read or write TKE file (en) in restart file 211 !! 212 !! ** Method : use of IOM library 213 !! if the restart does not contain TKE, en is either 214 !! set to rn_emin or recomputed 215 !!---------------------------------------------------------------------- 216 INTEGER , INTENT(in) :: kt ! ocean time-step 217 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 218 ! 219 INTEGER :: jit, jk ! dummy loop indices 220 INTEGER :: id1, id2 ! local integers 221 !!---------------------------------------------------------------------- 222 ! 223 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 224 ! ! --------------- 225 ! !* Read the restart file 226 IF( ln_rstart ) THEN 227 id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 228 id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 229 ! 230 IF( MIN( id1, id2 ) > 0 ) THEN ! restart exists => read it 231 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 232 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 233 ENDIF 234 ENDIF 235 ! !* otherwise Kz already set to the background value in zdf_phy_init 236 ! 237 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 238 ! ! ------------------- 239 IF(lwp) WRITE(numout,*) '---- ric-rst ----' 240 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 241 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 242 ! 243 ENDIF 244 ! 245 END SUBROUTINE ric_rst 213 246 214 247 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
