Changeset 11615 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
- Timestamp:
- 2019-09-30T15:18:21+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
r11329 r11615 2 2 !!====================================================================== 3 3 !! *** MODULE sbcblk_algo_coare3p0 *** 4 !! 5 !! After Fairall et al, 2003 4 6 !! Computes: 5 7 !! * bulk transfer coefficients C_D, C_E and C_H … … 7 9 !! * the effective bulk wind speed at 10m U_blk 8 10 !! => all these are used in bulk formulas in sbcblk.F90 9 !!10 !! Using the bulk formulation/param. of COARE v3, Fairall et al., 200311 11 !! 12 12 !! Routine turb_coare3p0 maintained and developed in AeroBulk … … 38 38 USE sbc_oce ! Surface boundary condition: ocean fields 39 39 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 USE sbcblk_skin 40 USE sbcblk_skin_coare3p6 ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 41 41 42 42 IMPLICIT NONE 43 43 PRIVATE 44 44 45 PUBLIC :: TURB_COARE3P0 ! called by sbcblk.F9045 PUBLIC :: COARE3P0_INIT, TURB_COARE3P0 46 46 47 47 ! !! COARE own values for given constants: … … 54 54 CONTAINS 55 55 56 SUBROUTINE turb_coare3p0( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 57 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 56 57 SUBROUTINE coare3p0_init(l_use_cs, l_use_wl) 58 !!--------------------------------------------------------------------- 59 !! *** FUNCTION coare3p0_init *** 60 !! 61 !! INPUT : 62 !! ------- 63 !! * l_use_cs : use the cool-skin parameterization 64 !! * l_use_wl : use the warm-layer parameterization 65 !!--------------------------------------------------------------------- 66 LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization 67 LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization 68 INTEGER :: ierr 69 !!--------------------------------------------------------------------- 70 IF ( l_use_wl ) THEN 71 ierr = 0 72 PRINT *, ' *** coare3p0_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj 73 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr ) 74 !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!' 75 Tau_ac(:,:) = 0._wp 76 Qnt_ac(:,:) = 0._wp 77 H_wl(:,:) = H_wl_max 78 PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!' 79 END IF 80 !! 81 IF ( l_use_cs ) THEN 82 ierr = 0 83 PRINT *, ' *** coare3p0_init: CS => allocating delta_vl :', jpi,jpj 84 ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr ) 85 !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of delta_vl and Qnt_ac failed!' 86 delta_vl(:,:) = 0.001_wp ! First guess of zdelta [m] 87 PRINT *, ' *** delta_vl allocated!' 88 END IF 89 END SUBROUTINE coare3p0_init 90 91 92 93 SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 94 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 58 95 & Cdn, Chn, Cen, & 59 & Qsw, rad_lw, slp, &60 & Tsk_b )96 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 97 & isecday_utc, plong, pdT_wl, Hwl ) ! optionals for warm-layer only 61 98 !!---------------------------------------------------------------------- 62 99 !! *** ROUTINE turb_coare3p0 *** … … 74 111 !! INPUT : 75 112 !! ------- 113 !! * kt : current time step (starts at 1) 76 114 !! * zt : height for temperature and spec. hum. of air [m] 77 115 !! * zu : height for wind speed (usually 10m) [m] 78 !! * U_zu : scalar wind speed at zu [m/s]79 116 !! * t_zt : potential air temperature at zt [K] 80 117 !! * q_zt : specific humidity of air at zt [kg/kg] 118 !! * U_zu : scalar wind speed at zu [m/s] 119 !! * l_use_cs : use the cool-skin parameterization 120 !! * l_use_wl : use the warm-layer parameterization 81 121 !! 82 122 !! INPUT/OUTPUT: … … 87 127 !! 88 128 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 89 !! -> doesn't need to be given a value if skin temp computed (in case l_use_ skin=True)90 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_ skin=False)91 !! 92 !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):129 !! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 130 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 131 !! 132 !! OPTIONAL INPUT: 93 133 !! --------------- 94 134 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 95 135 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 96 136 !! * slp : sea-level pressure [Pa] 97 !! * Tsk_b : estimate of skin temperature at previous time-step [K] 137 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 138 !! * isecday_utc: 139 !! * plong : longitude array [deg.E] 140 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 141 !! * Hwl : depth of warm layer [m] 98 142 !! 99 143 !! OUTPUT : … … 109 153 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 110 154 !!---------------------------------------------------------------------------------- 155 INTEGER, INTENT(in ) :: kt ! current time step 111 156 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 112 157 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] … … 116 161 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 117 162 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 163 LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization 164 LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization 118 165 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 119 166 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) … … 127 174 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 128 175 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 129 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Tsk_b ! [Pa] 176 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 177 ! 178 INTEGER, INTENT(in ), OPTIONAL :: isecday_utc ! current UTC time, counted in second since 00h of the current day 179 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: plong ! [deg.E] 180 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 181 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: Hwl ! [m] 130 182 ! 131 183 INTEGER :: j_itt … … 141 193 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt 142 194 ! 143 ! Cool skin: 144 LOGICAL :: l_use_skin = .FALSE. 145 REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST 195 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 196 & zsst, & ! to back up the initial bulk SST 197 & pdTc, & ! SST increment "dT" for cool-skin correction [K] 198 & pdTw, & ! SST increment "dT" for warm layer correction [K] 199 & zHwl ! depth of warm-layer [m] 200 201 ! 202 LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. 203 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0.F90' 204 CHARACTER(len=128) :: cf_tmp 146 205 !!---------------------------------------------------------------------------------- 147 206 148 ! Cool skin ? 149 IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 150 l_use_skin = .TRUE. 151 END IF 152 IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 207 IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays 208 153 209 154 210 l_zt_equal_zu = .FALSE. 155 211 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 156 157 212 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) 158 213 159 !! Initialization for cool skin: 160 zsst = T_s ! save the bulk SST 161 IF( l_use_skin ) THEN 162 ! First guess for skin temperature: 163 IF( PRESENT(Tsk_b) ) THEN 164 T_s = Tsk_b 165 ELSE 166 T_s = T_s - 0.25 ! sst - 0.25 167 END IF 168 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 169 END IF 214 !! Initializations for cool skin and warm layer: 215 IF ( l_use_cs ) THEN 216 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 217 PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!' 218 STOP 219 END IF 220 ALLOCATE ( pdTc(jpi,jpj) ) 221 pdTc(:,:) = -0.25_wp ! First guess of skin correction 222 END IF 223 224 IF ( l_use_wl ) THEN 225 IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN 226 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!' 227 STOP 228 END IF 229 ALLOCATE ( pdTw(jpi,jpj) ) 230 IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) ) 231 END IF 232 233 IF ( l_use_cs .OR. l_use_wl ) THEN 234 ALLOCATE ( zsst(jpi,jpj) ) 235 zsst = T_s ! backing up the bulk SST 236 IF( l_use_cs ) T_s = T_s - 0.25 ! First guess of correction 237 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 238 END IF 239 170 240 171 241 !! First guess of temperature and humidity at height zu: … … 190 260 z0t = MIN(ABS(z0t), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 191 261 192 ztmp2 = vkarmn/ztmp0 193 Cd = ztmp2*ztmp2 ! first guess of Cd 262 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 194 263 195 264 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd … … 234 303 ztmp1 = u_star*u_star ! u*^2 235 304 236 !! Update wind at zu taking into acount convection-related wind gustiness: 237 ! Ug = Beta*w* (Beta = 1.25, Fairall et al. 2003, Eq.8): 238 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2 239 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 240 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 305 !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8): 306 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 307 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 308 U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 241 309 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 242 310 … … 251 319 !! Adjustment the wind at 10m (not needed in the current algo form): 252 320 !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 253 321 254 322 !! Roughness lengthes z0, z0t (z0q = z0t) : 255 323 ztmp2 = u_star/vkarmn*LOG(10./z0) ! Neutral wind speed at 10m 256 324 z0 = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) 257 ztmp1 = z0*u_star/znu_a ! Re_r: roughness Reynolds number258 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 **(-0.6_wp)) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!325 ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 326 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!! 259 327 260 328 !! Turbulent scales at zu : … … 267 335 268 336 IF( .NOT. l_zt_equal_zu ) THEN 269 ! What's need to be done if zt /= zu 270 !! Re-updating temperature and humidity at zu : 271 ztmp2 = ztmp0 - psi_h_coare(zeta_t) 272 ztmp1 = log(zt/zu) + ztmp2 337 !! Re-updating temperature and humidity at zu if zt /= zu : 338 ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t) 273 339 t_zu = t_zt - t_star/vkarmn*ztmp1 274 340 q_zu = q_zt - q_star/vkarmn*ztmp1 275 341 END IF 276 342 277 !! SKIN related part 278 !! ----------------- 279 IF( l_use_skin ) THEN 280 !! compute transfer coefficients at zu : lolo: verifier... 281 ztmp0 = u_star/U_blk 282 Ch = ztmp0*t_star/dt_zu 283 Ce = ztmp0*q_star/dq_zu 284 ! Non-Solar heat flux to the ocean: 285 ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp) ! rho*U10 286 ztmp2 = T_s*T_s 287 ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux 288 & + emiss_w*(rad_lw - stefan*ztmp2*ztmp2) ! Net longwave flux 289 !! => "ztmp1" is the net non-solar surface heat flux ! 290 !! Updating the values of the skin temperature T_s and q_s : 291 CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!) 292 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! 200 -> just to avoid numerics problem on masked regions if silly values are given 293 END IF 294 295 IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 343 344 IF( l_use_cs ) THEN 345 !! Cool-skin contribution 346 347 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 348 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 349 350 CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2, pdTc ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 351 352 T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1) 353 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1) 354 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 355 356 END IF 357 358 IF( l_use_wl ) THEN 359 !! Warm-layer contribution 360 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, & 361 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 362 !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 363 IF (PRESENT(Hwl)) THEN 364 CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt), pdTw, Hwl=zHwl ) 365 ELSE 366 CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt), pdTw ) 367 END IF 368 !! Updating T_s and q_s !!! 369 T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1) 370 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1) 371 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 372 !LOLO: 373 !IF ( j_itt == nb_itt) THEN 374 ! WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 375 ! CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 376 ! WRITE(cf_tmp, '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 377 ! CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 378 !END IF 379 !LOLO. 380 END IF 381 382 383 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 296 384 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 297 385 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) … … 312 400 313 401 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 402 403 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 404 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 405 406 IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 407 IF ( l_use_cs ) DEALLOCATE ( pdTc ) 408 IF ( l_use_wl ) THEN 409 DEALLOCATE ( pdTw ) 410 IF (PRESENT(Hwl)) DEALLOCATE ( zHwl ) 411 END IF 314 412 315 413 END SUBROUTINE turb_coare3p0
Note: See TracChangeset
for help on using the changeset viewer.
