Changeset 15084
- Timestamp:
- 2021-07-05T21:48:42+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/ISF
- Files:
-
- 3 edited
-
isfcav.F90 (modified) (8 diffs)
-
isfcavgam.F90 (modified) (11 diffs)
-
isfpar.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ISF/isfcav.F90
r15082 r15084 22 22 USE isfdiags , ONLY: isf_diags_flx ! ice shelf diags subroutine 23 23 ! 24 USE oce , ONLY: ts ! oceantracers24 USE oce , ONLY: ts, uu, vv, rn2 ! ocean dynamics and tracers 25 25 USE par_oce , ONLY: jpi,jpj ! ocean space and time domain 26 26 USE phycst , ONLY: grav,rho0,rho0_rcp,r1_rho0_rcp ! physical constants … … 31 31 USE fldread ! read input field at current time step 32 32 USE lbclnk ! lbclnk 33 USE lib_mpp ! MPP library 33 34 34 35 IMPLICIT NONE … … 38 39 PUBLIC isf_cav, isf_cav_init ! routine called in isfmlt 39 40 41 !! * Substitutions 42 # include "do_loop_substitute.h90" 40 43 !!---------------------------------------------------------------------- 41 44 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 71 74 !!--------------------------------------------------------------------- 72 75 LOGICAL :: lit 73 INTEGER :: nit 76 INTEGER :: nit, ji, jj, ikt 74 77 REAL(wp) :: zerr 78 REAL(wp) :: zcoef, zdku, zdkv 75 79 REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh ! heat fluxes 76 REAL(wp), DIMENSION(jpi,jpj) :: zqh_b 80 REAL(wp), DIMENSION(jpi,jpj) :: zqh_b, zRc ! 77 81 REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas ! exchange coeficient 78 82 REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl ! temp. and sal. in top boundary layer … … 91 95 zqoce(:,:) = -pqfwf(:,:) * rLfusisf ! 92 96 zqh_b(:,:) = ptsc(:,:,jp_tem) * rho0_rcp ! last time step total heat fluxes (to speed up convergence) 97 98 DO_2D( 0, 0, 0, 0 ) 99 ikt = mikt(ji,jj) 100 ! compute Rc number (as done in zdfric.F90) 101 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 102 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 103 ! ! shear of horizontal velocity 104 zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) & 105 & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) ) 106 zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) & 107 & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) ) 108 ! ! richardson number (minimum value set to zero) 109 zRc(ji,jj) = MAX(rn2(ji,jj,ikt+1), 1.e-20_wp) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp ) 110 END_2D 111 CALL lbc_lnk( 'isfmlt', zRc, 'T', 1._wp ) 93 112 ENDIF 94 113 ! … … 100 119 ! useless if melt specified 101 120 IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN 102 CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce , pqfwf, &121 CALL isfcav_gammats( Kmm, zttbl, zstbl, zqoce , pqfwf, zRc, & 103 122 & zgammat, zgammas ) 104 123 END IF … … 115 134 CASE ( 'vel_stab' ) 116 135 ! compute error between 2 iterations 117 zerr = MAXVAL(ABS(zqhc(:,:)+zqoce(:,:) - zqh_b(:,:))) 136 zerr = 0._wp 137 DO_2D( 0, 0, 0, 0 ) 138 zerr = MAX( zerr, ABS(zqhc(ji,jj)+zqoce(ji,jj) - zqh_b(ji,jj)) ) 139 END_2D 140 CALL mpp_max( 'isfcav', zerr ) ! max over the global domain 118 141 ! 119 142 ! define if iteration needed … … 130 153 END DO 131 154 ! 132 ! compute heat and water flux ( > 0 from isf to oce) 133 pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:) 134 zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:) 135 zqhc (:,:) = zqhc(:,:) * mskisf_cav(:,:) 136 ! 137 ! compute heat content flux ( > 0 from isf to oce) 138 zqlat(:,:) = - pqfwf(:,:) * rLfusisf ! 2d latent heat flux (W/m2) 139 ! 140 ! total heat flux ( > 0 from isf to oce) 141 zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) ) 142 ! 143 ! lbclnk on melt 144 CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp) 155 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 156 ! compute heat and water flux ( > 0 from isf to oce) 157 pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_cav(ji,jj) 158 zqoce(ji,jj) = zqoce(ji,jj) * mskisf_cav(ji,jj) 159 zqhc (ji,jj) = zqhc(ji,jj) * mskisf_cav(ji,jj) 160 ! 161 ! compute heat content flux ( > 0 from isf to oce) 162 zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2) 163 ! 164 ! total heat flux ( > 0 from isf to oce) 165 zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) ) 166 ! 167 ! set temperature content 168 ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp 169 END_2D 145 170 ! 146 171 ! output fluxes 147 172 CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc) 148 !149 ! set temperature content150 ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp151 173 ! 152 174 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) -
NEMO/trunk/src/OCE/ISF/isfcavgam.F90
r15082 r15084 14 14 USE isftbl , ONLY: isf_tbl 15 15 16 USE oce , ONLY: uu, vv , rn2 ! ocean dynamics and tracers16 USE oce , ONLY: uu, vv ! ocean dynamics 17 17 USE phycst , ONLY: grav, vkarmn ! physical constant 18 18 USE eosbn2 , ONLY: eos_rab ! equation of state … … 44 44 !!----------------------------------------------------------------------------------------------------- 45 45 ! 46 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, p gt, pgs )46 SUBROUTINE isfcav_gammats( Kmm, pttbl, pstbl, pqoce, pqfwf, pRc, pgt, pgs ) 47 47 !!---------------------------------------------------------------------- 48 48 !! ** Purpose : compute the coefficient echange for heat and fwf flux … … 57 57 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqoce, pqfwf ! isf heat and fwf 58 58 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer 59 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number 59 60 !!--------------------------------------------------------------------- 60 61 REAL(wp), DIMENSION(jpi,jpj) :: zutbl, zvtbl ! top boundary layer velocity … … 94 95 pgs(:,:) = rn_gammas0 95 96 CASE ( 'vel' ) ! gamma is proportional to u* 96 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs )97 CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs ) 97 98 CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u* 98 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, p gt, pgs )99 CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs ) 99 100 CASE DEFAULT 100 101 CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)') … … 153 154 END SUBROUTINE gammats_vel 154 155 155 SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, & ! <<== in156 & pgt , pgs ) ! ==>> out gammats [m/s]156 SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, pRc, & ! <<== in 157 & pgt , pgs ) ! ==>> out gammats [m/s] 157 158 !!---------------------------------------------------------------------- 158 159 !! ** Purpose : compute the coefficient echange coefficient … … 171 172 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: putbl, pvtbl ! velocity in the losch top boundary layer 172 173 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! tracer in the losch top boundary layer 174 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pRc ! Richardson number 173 175 !!--------------------------------------------------------------------- 174 176 INTEGER :: ji, jj ! loop index 175 177 INTEGER :: ikt ! local integer 176 178 REAL(wp) :: zdku, zdkv ! U, V shear 177 REAL(wp) :: zPr, zSc , zRc ! Prandtl, Scmidth and Richardsonnumber179 REAL(wp) :: zPr, zSc ! Prandtl and Scmidth number 178 180 REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point 179 181 REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness … … 190 192 !!--------------------------------------------------------------------- 191 193 ! 192 ! compute ustar193 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )194 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) )195 END_2D196 !197 ! output ustar198 CALL iom_put('isfustar',zustar(:,:))199 !200 194 ! compute Pr and Sc number (eq ??) 201 195 zPr = 13.8_wp … … 207 201 ! 208 202 ! compute gamma 209 DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 203 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 204 210 205 ikt = mikt(ji,jj) 206 207 ! compute ustar 208 zustar(ji,jj) = SQRT( pCd(ji,jj) * ( putbl(ji,jj) * putbl(ji,jj) + pvtbl(ji,jj) * pvtbl(ji,jj) + pke2 ) ) 211 209 212 210 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think … … 214 212 pgs(ji,jj) = rn_gammas0 215 213 ELSE 216 ! compute Rc number (as done in zdfric.F90)217 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation218 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)219 ! ! shear of horizontal velocity220 zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) &221 & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) )222 zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) &223 & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) )224 ! ! richardson number (minimum value set to zero)225 zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )226 227 214 ! compute bouyancy 228 215 zts(jp_tem) = pttbl(ji,jj) … … 244 231 ! 245 232 ! compute eta* (stability parameter) (Eq ??) 246 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 233 zetastar = 1._wp / ( SQRT(1._wp + MAX( 0._wp, zxsiN * zustar(ji,jj) & 234 & / MAX( 1.e-20, ABS(ff_t(ji,jj)) * zmols * pRc(ji,jj) ) ))) 247 235 ! 248 236 ! compute the sublayer thickness (Eq ??) 249 zhnu = 5 * znu / zustar(ji,jj)237 zhnu = 5 * znu / MAX( 1.e-20, zustar(ji,jj) ) 250 238 ! 251 239 ! compute gamma turb (Eq ??) 252 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &240 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / MAX( 1.e-10, ABS(ff_t(ji,jj)) * zhnu )) & 253 241 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 254 242 ! … … 258 246 END IF 259 247 END_2D 248 ! output ustar 249 CALL iom_put('isfustar',zustar(:,:)) 260 250 261 251 END SUBROUTINE gammats_vel_stab -
NEMO/trunk/src/OCE/ISF/isfpar.F90
r15004 r15084 30 30 USE iom ! I/O library 31 31 USE fldread ! read input field at current time step 32 USE lbclnk ! lbc_lnk33 32 34 33 IMPLICIT NONE … … 37 36 PUBLIC isf_par, isf_par_init 38 37 38 !! * Substitutions 39 # include "do_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 41 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 67 68 INTEGER, INTENT(in) :: Kmm ! ocean time level index 68 69 !!--------------------------------------------------------------------- 70 INTEGER :: ji, jj 69 71 REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh 70 72 !!--------------------------------------------------------------------- … … 73 75 CALL isfpar_mlt( kt, Kmm, zqhc, zqoce, pqfwf ) 74 76 ! 75 ! compute heat and water flux (from isf to oce) 76 pqfwf(:,:) = pqfwf(:,:) * mskisf_par(:,:) 77 zqoce(:,:) = zqoce(:,:) * mskisf_par(:,:) 78 zqhc (:,:) = zqhc(:,:) * mskisf_par(:,:) 79 ! 80 ! compute latent heat flux (from isf to oce) 81 zqlat(:,:) = - pqfwf(:,:) * rLfusisf ! 2d latent heat flux (W/m2) 82 ! 83 ! total heat flux (from isf to oce) 84 zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) ) 85 ! 86 ! lbclnk on melt and heat fluxes 87 CALL lbc_lnk( 'isfmlt', zqh, 'T', 1.0_wp, pqfwf, 'T', 1.0_wp) 77 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 78 ! compute heat and water flux (from isf to oce) 79 pqfwf(ji,jj) = pqfwf(ji,jj) * mskisf_par(ji,jj) 80 zqoce(ji,jj) = zqoce(ji,jj) * mskisf_par(ji,jj) 81 zqhc (ji,jj) = zqhc(ji,jj) * mskisf_par(ji,jj) 82 ! 83 ! compute latent heat flux (from isf to oce) 84 zqlat(ji,jj) = - pqfwf(ji,jj) * rLfusisf ! 2d latent heat flux (W/m2) 85 ! 86 ! total heat flux (from isf to oce) 87 zqh(ji,jj) = ( zqhc (ji,jj) + zqoce(ji,jj) ) 88 ! 89 ! set temperature content 90 ptsc(ji,jj,jp_tem) = zqh(ji,jj) * r1_rho0_rcp 91 END_2D 88 92 ! 89 93 ! output fluxes 90 94 CALL isf_diags_flx( Kmm, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc) 91 !92 ! set temperature content93 ptsc(:,:,jp_tem) = zqh(:,:) * r1_rho0_rcp94 95 ! 95 96 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
Note: See TracChangeset
for help on using the changeset viewer.
