- Timestamp:
- 2017-05-30T10:13:14+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r7991 r8093 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 27 USE eosbn2 , ONLY : neos28 27 ! 29 28 USE in_out_manager ! I/O manager 30 USE lbclnk ! ocean lateral boundary condition (or mpp link) 31 USE lib_mpp ! MPP library 29 USE iom ! I/O manager library 32 30 USE timing ! Timing 33 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) … … 37 35 PRIVATE 38 36 39 PUBLIC zdf_ric ! called by step.F90 40 PUBLIC zdf_ric_init ! called by opa.F90 37 PUBLIC zdf_ric ! called by zdfphy.F90 38 PUBLIC ric_rst ! called by zdfphy.F90 39 PUBLIC zdf_ric_init ! called by nemogcm.F90 41 40 42 41 ! !!* Namelist namzdf_ric : Richardson number dependent Kz * … … 54 53 # include "vectopt_loop_substitute.h90" 55 54 !!---------------------------------------------------------------------- 56 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)55 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 57 56 !! $Id$ 58 57 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 59 58 !!---------------------------------------------------------------------- 60 59 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 60 162 61 SUBROUTINE zdf_ric_init … … 205 104 ENDIF 206 105 ! 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) 106 CALL ric_rst( nit000, 'READ' ) !* read or initialize all required files 107 ! 108 END SUBROUTINE zdf_ric_init 109 110 111 SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt ) 112 !!---------------------------------------------------------------------- 113 !! *** ROUTINE zdfric *** 114 !! 115 !! ** Purpose : Compute the before eddy viscosity and diffusivity as 116 !! a function of the local richardson number. 117 !! 118 !! ** Method : Local richardson number dependent formulation of the 119 !! vertical eddy viscosity and diffusivity coefficients. 120 !! The eddy coefficients are given by: 121 !! avm = avm0 + avmb 122 !! avt = avm0 / (1 + rn_alp*ri) 123 !! with ri = N^2 / dz(u)**2 124 !! = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ] 125 !! avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 126 !! where ri is the before local Richardson number, 127 !! rn_avmri is the maximum value reaches by avm and avt 128 !! and rn_alp, nn_ric are adjustable parameters. 129 !! Typical values : rn_alp=5. and nn_ric=2. 130 !! 131 !! As second step compute Ekman depth from wind stress forcing 132 !! and apply namelist provided vertical coeff within this depth. 133 !! The Ekman depth is: 134 !! Ustar = SQRT(Taum/rho0) 135 !! ekd= rn_ekmfc * Ustar / f0 136 !! Large et al. (1994, eq.24) suggest rn_ekmfc=0.7; however, the derivation 137 !! of the above equation indicates the value is somewhat arbitrary; therefore 138 !! we allow the freedom to increase or decrease this value, if the 139 !! Ekman depth estimate appears too shallow or too deep, respectively. 140 !! Ekd is then limited by rn_mldmin and rn_mldmax provided in the 141 !! namelist 142 !! N.B. the mask are required for implicit scheme, and surface 143 !! and bottom value already set in zdfphy.F90 144 !! 145 !! ** Action : avm, avt mixing coeff (inner domain values only) 146 !! 147 !! References : Pacanowski & Philander 1981, JPO, 1441-1451. 148 !! PFJ Lermusiaux 2001. 149 !!---------------------------------------------------------------------- 150 INTEGER , INTENT(in ) :: kt ! ocean time-step 151 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdept ! depth of t-point [m] 152 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 153 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 154 !! 155 INTEGER :: ji, jj, jk ! dummy loop indices 156 REAL(wp) :: zcfRi, zav, zustar, zhek ! local scalars 157 REAL(wp), DIMENSION(jpi,jpj) :: zh_ekm ! 2D workspace 158 !!---------------------------------------------------------------------- 159 ! 160 IF( nn_timing == 1 ) CALL timing_start('zdf_ric') 161 ! 162 ! !== avm and avt = F(Richardson number) ==! 163 DO jk = 2, jpkm1 164 DO jj = 1, jpjm1 165 DO ji = 1, jpim1 ! coefficient = F(richardson number) (avm-weighted Ri) 166 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 ) ) ) 167 zav = rn_avmri * zcfRi**nn_ric 168 ! ! avm and avt coefficients 169 p_avm(ji,jj,jk) = MAX( zav , avmb(jk) ) * wmask(ji,jj,jk) 170 p_avt(ji,jj,jk) = MAX( zav * zcfRi , avtb(jk) ) * wmask(ji,jj,jk) 171 END DO 172 END DO 210 173 END DO 211 174 ! 212 END SUBROUTINE zdf_ric_init 175 !!gm BUG <<<<==== This param can't work at low latitude 176 !!gm it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! ) 177 ! 178 IF( ln_mldw ) THEN !== set a minimum value in the Ekman layer ==! 179 ! 180 DO jj = 2, jpjm1 !* Ekman depth 181 DO ji = 2, jpim1 182 zustar = SQRT( taum(ji,jj) * r1_rau0 ) 183 zhek = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall ) ! Ekman depth 184 zh_ekm(ji,jj) = MAX( rn_mldmin , MIN( zhek , rn_mldmax ) ) ! set allowed range 185 END DO 186 END DO 187 DO jk = 2, jpkm1 !* minimum mixing coeff. within the Ekman layer 188 DO jj = 2, jpjm1 189 DO ji = 2, jpim1 190 IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN 191 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) 192 p_avt(ji,jj,jk) = MAX( p_avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk) 193 ENDIF 194 END DO 195 END DO 196 END DO 197 ENDIF 198 ! 199 IF( nn_timing == 1 ) CALL timing_stop('zdf_ric') 200 ! 201 END SUBROUTINE zdf_ric 202 203 204 SUBROUTINE ric_rst( kt, cdrw ) 205 !!--------------------------------------------------------------------- 206 !! *** ROUTINE ric_rst *** 207 !! 208 !! ** Purpose : Read or write TKE file (en) in restart file 209 !! 210 !! ** Method : use of IOM library 211 !! if the restart does not contain TKE, en is either 212 !! set to rn_emin or recomputed 213 !!---------------------------------------------------------------------- 214 INTEGER , INTENT(in) :: kt ! ocean time-step 215 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 216 ! 217 INTEGER :: jit, jk ! dummy loop indices 218 INTEGER :: id1, id2 ! local integers 219 !!---------------------------------------------------------------------- 220 ! 221 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 222 ! ! --------------- 223 ! !* Read the restart file 224 IF( ln_rstart ) THEN 225 id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 226 id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 227 ! 228 IF( MIN( id1, id2 ) > 0 ) THEN ! restart exists => read it 229 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 230 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 231 ENDIF 232 ENDIF 233 ! !* otherwise Kz already set to the background value in zdf_phy_init 234 ! 235 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 236 ! ! ------------------- 237 IF(lwp) WRITE(numout,*) '---- ric-rst ----' 238 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 239 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 240 ! 241 ENDIF 242 ! 243 END SUBROUTINE ric_rst 213 244 214 245 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
