- Timestamp:
- 2019-10-23T16:04:12+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90
r11615 r11772 82 82 REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3] 83 83 REAL(wp), PARAMETER, PUBLIC :: rho0_w = 1025._wp !: Density of sea-water (ECMWF->1025) [kg/m^3] 84 REAL(wp), PARAMETER, PUBLIC :: roadrw = rho0_a/rho0_w !: Density ratio 84 REAL(wp), PARAMETER, PUBLIC :: radrw = rho0_a/rho0_w !: Density ratio 85 REAL(wp), PARAMETER, PUBLIC :: sq_radrw = SQRT(rho0_a/rho0_w) 85 86 REAL(wp), PARAMETER, PUBLIC :: rCp0_w = 4190._wp !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg] 86 87 REAL(wp), PARAMETER, PUBLIC :: rnu0_w = 1.e-6_wp !: kinetic viscosity of water [m^2/s] -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r11672 r11772 542 542 CASE( np_ECMWF ) 543 543 IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 544 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF544 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF 545 545 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 546 546 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) … … 585 585 CASE( np_ECMWF ) 586 586 IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' !LBrm 587 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF587 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & ! ECMWF 588 588 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 589 589 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
r11672 r11772 49 49 REAL(wp), PARAMETER :: Beta0 = 1.25_wp ! gustiness parameter 50 50 51 INTEGER , PARAMETER :: nb_itt = 5! number of itterations51 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 52 52 53 53 !!---------------------------------------------------------------------- … … 70 70 IF ( l_use_wl ) THEN 71 71 ierr = 0 72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H _wl(jpi,jpj), STAT=ierr )72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), Hz_wl(jpi,jpj), dT_wl(jpi,jpj), STAT=ierr ) 73 73 !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!' 74 74 Tau_ac(:,:) = 0._wp 75 Qnt_ac(:,:) = 0._wp 76 H_wl(:,:) = H_wl_max 75 Qnt_ac(:,:) = 0._wp 76 Hz_wl(:,:) = Hwl_max 77 dT_wl(:,:) = 0._wp 77 78 END IF 78 79 !! 79 80 IF ( l_use_cs ) THEN 80 81 ierr = 0 81 ALLOCATE ( d elta_vl(jpi,jpj), STAT=ierr )82 !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of d elta_vland Qnt_ac failed!'83 d elta_vl(:,:) = 0.001_wp ! First guess of zdelta [m]82 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 83 !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of dT_cs and Qnt_ac failed!' 84 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 84 85 END IF 85 86 END SUBROUTINE coare3p0_init … … 91 92 & Cdn, Chn, Cen, & 92 93 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 93 & pdT_wl, Hwl ) ! optionals for warm-layer only94 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 94 95 !!---------------------------------------------------------------------- 95 96 !! *** ROUTINE turb_coare3p0 *** … … 133 134 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 134 135 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 135 !! * Hwl : depth of warm layer[m]136 !! * pHz_wl : thickness of warm-layer [m] 136 137 !! 137 138 !! OUTPUT : … … 171 172 ! 172 173 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 173 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: Hwl! [m]174 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 174 175 ! 175 176 INTEGER :: j_itt … … 190 191 & pdTw, & ! SST increment "dT" for warm layer correction [K] 191 192 & zHwl ! depth of warm-layer [m] 193 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0' 192 194 !!---------------------------------------------------------------------------------- 193 195 194 IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays196 IF ( kt == nit000 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) 195 197 196 198 l_zt_equal_zu = .FALSE. … … 199 201 200 202 !! Initializations for cool skin and warm layer: 203 IF ( l_use_cs ) THEN 204 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 205 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 206 END IF 207 END IF 208 209 IF ( l_use_wl ) THEN 210 IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) )) THEN 211 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 212 END IF 213 END IF 214 201 215 IF ( l_use_cs .OR. l_use_wl ) THEN 202 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) &203 & CALL ctl_stop( 'turb_coare3p0 => provide Qsw, rad_lw & slp to ', 'use cool-skin and/or warm-layer param' )204 216 ALLOCATE ( zsst(jpi,jpj) ) 205 217 zsst = T_s ! backing up the bulk SST 206 IF( l_use_cs ) T_s = T_s - 0.25 ! First guess of correction218 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 207 219 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 208 END IF209 IF ( l_use_cs ) THEN210 ALLOCATE ( pdTc(jpi,jpj) )211 pdTc(:,:) = -0.25_wp ! First guess of skin correction212 END IF213 IF ( l_use_wl ) THEN214 ALLOCATE ( pdTw(jpi,jpj) )215 IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) )216 220 END IF 217 221 … … 323 327 !! Cool-skin contribution 324 328 325 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &329 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 326 330 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 327 331 328 CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 , pdTc) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2329 330 T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)331 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1)332 CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 333 334 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 335 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 332 336 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 333 337 … … 336 340 IF( l_use_wl ) THEN 337 341 !! Warm-layer contribution 338 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &342 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 339 343 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 340 344 !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 341 IF (PRESENT(Hwl)) THEN 342 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt), pdTw, Hwl=zHwl ) 343 ELSE 344 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt), pdTw ) 345 END IF 345 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt) ) 346 346 347 !! Updating T_s and q_s !!! 347 T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1)348 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1)348 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 349 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 349 350 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 350 !LOLO: 351 !IF ( j_itt == nb_itt) THEN 352 ! WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 353 ! CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 354 ! WRITE(cf_tmp, '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 355 ! CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 356 !END IF 357 !LOLO. 358 END IF 359 351 END IF 360 352 361 353 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN … … 379 371 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 380 372 381 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 382 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 373 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 374 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 375 IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 383 376 384 377 IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 385 IF ( l_use_cs ) DEALLOCATE ( pdTc )386 IF ( l_use_wl ) THEN387 DEALLOCATE ( pdTw )388 IF (PRESENT(Hwl)) DEALLOCATE ( zHwl )389 END IF390 378 391 379 END SUBROUTINE turb_coare3p0 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
r11672 r11772 45 45 PUBLIC :: COARE3P6_INIT, TURB_COARE3P6 46 46 47 ! !! COARE own values for given constants:48 REAL(wp), PARAMETER :: zi0= 600._wp ! scale height of the atmospheric boundary layer...49 REAL(wp), PARAMETER :: Beta0 = 1.2_wp! gustiness parameter50 51 INTEGER , PARAMETER :: nb_itt = 5! number of itterations47 !! COARE own values for given constants: 48 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 49 REAL(wp), PARAMETER :: Beta0 = 1.2_wp ! gustiness parameter 50 51 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 52 52 53 53 !!---------------------------------------------------------------------- … … 70 70 IF ( l_use_wl ) THEN 71 71 ierr = 0 72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H _wl(jpi,jpj), STAT=ierr )73 !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac and Qnt_acfailed!'72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), Hz_wl(jpi,jpj), dT_wl(jpi,jpj), STAT=ierr ) 73 !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of Tau_ac, Qnt_ac, dT_wl & Hz_wl failed!' 74 74 Tau_ac(:,:) = 0._wp 75 Qnt_ac(:,:) = 0._wp 76 H_wl(:,:) = H_wl_max 75 Qnt_ac(:,:) = 0._wp 76 dT_wl(:,:) = 0._wp 77 Hz_wl(:,:) = Hwl_max 77 78 END IF 78 79 !! 79 80 IF ( l_use_cs ) THEN 80 81 ierr = 0 81 ALLOCATE ( d elta_vl(jpi,jpj), STAT=ierr )82 !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of d elta_vl and Qnt_acfailed!'83 d elta_vl(:,:) = 0.001_wp ! First guess of zdelta [m]82 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 83 !IF( ierr > 0 ) STOP ' COARE3P6_INIT => allocation of dT_cs failed!' 84 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 84 85 END IF 85 86 END SUBROUTINE coare3p6_init … … 91 92 & Cdn, Chn, Cen, & 92 93 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 93 & pdT_wl, Hwl ) ! optionals for warm-layer only94 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 94 95 !!---------------------------------------------------------------------- 95 96 !! *** ROUTINE turb_coare3p6 *** … … 131 132 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 132 133 !! * slp : sea-level pressure [Pa] 134 !! 135 !! OPTIONAL OUTPUT: 136 !! ---------------- 133 137 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 134 138 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 135 !! * Hwl : depth of warm layer[m]139 !! * pHz_wl : thickness of warm-layer [m] 136 140 !! 137 141 !! OUTPUT : … … 169 173 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 170 174 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 171 !172 175 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 173 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: Hwl! [m]176 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 174 177 ! 175 178 INTEGER :: j_itt … … 190 193 & pdTw, & ! SST increment "dT" for warm layer correction [K] 191 194 & zHwl ! depth of warm-layer [m] 195 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6' 192 196 !!---------------------------------------------------------------------------------- 193 197 194 IF ( kt == 1 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays198 IF ( kt == nit000 ) CALL COARE3P6_INIT(l_use_cs, l_use_wl) 195 199 196 200 l_zt_equal_zu = .FALSE. … … 199 203 200 204 !! Initializations for cool skin and warm layer: 205 IF ( l_use_cs ) THEN 206 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 207 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 208 END IF 209 END IF 210 211 IF ( l_use_wl ) THEN 212 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 213 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 214 END IF 215 END IF 216 201 217 IF ( l_use_cs .OR. l_use_wl ) THEN 202 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) &203 & CALL ctl_stop( 'turb_coare3p6 => provide Qsw, rad_lw & slp to ', 'use cool-skin and/or warm-layer param' )204 218 ALLOCATE ( zsst(jpi,jpj) ) 205 219 zsst = T_s ! backing up the bulk SST 206 IF( l_use_cs ) T_s = T_s - 0.25 ! First guess of correction220 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 207 221 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!! 208 END IF209 IF ( l_use_cs ) THEN210 ALLOCATE ( pdTc(jpi,jpj) )211 pdTc(:,:) = -0.25_wp ! First guess of skin correction212 END IF213 IF ( l_use_wl ) THEN214 ALLOCATE ( pdTw(jpi,jpj) )215 IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) )216 222 END IF 217 223 … … 323 329 !! Cool-skin contribution 324 330 325 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &331 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 326 332 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 327 333 328 CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 , pdTc) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2329 330 T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)331 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1)334 CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 ) ! ! Qnsol -> ztmp1 / Qlat -> ztmp2 335 336 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 337 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 332 338 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 333 339 … … 336 342 IF( l_use_wl ) THEN 337 343 !! Warm-layer contribution 338 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &344 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 339 345 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 340 346 !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 341 IF (PRESENT(Hwl)) THEN 342 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt), pdTw, Hwl=zHwl ) 343 ELSE 344 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt), pdTw ) 345 END IF 347 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt) ) 348 346 349 !! Updating T_s and q_s !!! 347 T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1)348 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1)350 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 351 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 349 352 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 350 !LOLO: 351 !IF ( j_itt == nb_itt) THEN 352 ! WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea 353 ! CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac') 354 ! WRITE(cf_tmp, '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea 355 ! CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw') 356 !END IF 357 !LOLO. 358 END IF 359 353 END IF 360 354 361 355 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN … … 379 373 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 380 374 381 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1) 382 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1) 375 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 376 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 377 IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 383 378 384 379 IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 385 IF ( l_use_cs ) DEALLOCATE ( pdTc )386 IF ( l_use_wl ) THEN387 DEALLOCATE ( pdTw )388 IF (PRESENT(Hwl)) DEALLOCATE ( zHwl )389 END IF390 380 391 381 END SUBROUTINE turb_coare3p6 … … 395 385 !!------------------------------------------------------------------- 396 386 !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m 397 !! 387 !! "wind speed dependent formulation" 398 388 !! (Eq. 13 in Edson et al., 2013) 399 389 !! … … 405 395 REAL(wp), PARAMETER :: charn0_max = 0.028 !: value above which the Charnock parameter levels off for winds > 18 m/s 406 396 !!------------------------------------------------------------------- 407 !408 397 alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp ) 409 ! 398 !! 410 399 END FUNCTION alfa_charn_3p6 400 401 FUNCTION alfa_charn_3p6_wave( pus, pwsh, pwps ) 402 !!------------------------------------------------------------------- 403 !! Computes the Charnock parameter as a function of wave information and u* 404 !! 405 !! (COARE 3.6, Fairall et al., 2018) 406 !! 407 !! Author: L. Brodeau, October 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 408 !!------------------------------------------------------------------- 409 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6_wave 410 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! friction velocity [m/s] 411 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwsh ! significant wave height [m] 412 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwps ! phase speed of dominant waves [m/s] 413 !!------------------------------------------------------------------- 414 alfa_charn_3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus) 415 !! 416 END FUNCTION alfa_charn_3p6_wave 417 411 418 412 419 FUNCTION psi_m_coare( pzeta ) -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r11631 r11772 44 44 PRIVATE 45 45 46 PUBLIC :: TURB_ECMWF ! called by sbcblk.F9046 PUBLIC :: ECMWF_INIT, TURB_ECMWF 47 47 48 48 ! !! ECMWF own values for given constants, taken form IFS documentation... … … 55 55 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 56 56 57 INTEGER , PARAMETER :: nb_itt = 5! number of itterations57 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 58 58 59 59 !!---------------------------------------------------------------------- 60 60 CONTAINS 61 61 62 SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 63 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 64 & Cdn, Chn, Cen, & 65 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 66 & pdT_wl ) ! optionals for warm-layer only 62 63 SUBROUTINE ecmwf_init(l_use_cs, l_use_wl) 64 !!--------------------------------------------------------------------- 65 !! *** FUNCTION ecmwf_init *** 66 !! 67 !! INPUT : 68 !! ------- 69 !! * l_use_cs : use the cool-skin parameterization 70 !! * l_use_wl : use the warm-layer parameterization 71 !!--------------------------------------------------------------------- 72 LOGICAL , INTENT(in) :: l_use_cs ! use the cool-skin parameterization 73 LOGICAL , INTENT(in) :: l_use_wl ! use the warm-layer parameterization 74 INTEGER :: ierr 75 !!--------------------------------------------------------------------- 76 IF ( l_use_wl ) THEN 77 ierr = 0 78 ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 79 !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' 80 dT_wl(:,:) = 0._wp 81 Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 82 END IF 83 !! 84 IF ( l_use_cs ) THEN 85 ierr = 0 86 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 87 !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_cs failed!' 88 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 89 END IF 90 END SUBROUTINE ecmwf_init 91 92 93 94 SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 95 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 96 & Cdn, Chn, Cen, & 97 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 98 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 67 99 !!---------------------------------------------------------------------- 68 100 !! *** ROUTINE turb_ecmwf *** … … 80 112 !! INPUT : 81 113 !! ------- 114 !! * kt : current time step (starts at 1) 82 115 !! * zt : height for temperature and spec. hum. of air [m] 83 116 !! * zu : height for wind speed (usually 10m) [m] … … 95 128 !! 96 129 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 97 !! -> doesn't need to be given a value if skin temp computed (in case l_use_ skin=True)98 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_ skin=False)130 !! -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True) 131 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False) 99 132 !! 100 133 !! OPTIONAL INPUT: … … 103 136 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 104 137 !! * slp : sea-level pressure [Pa] 138 !! 139 !! OPTIONAL OUTPUT: 140 !! ---------------- 105 141 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 106 142 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 143 !! * pHz_wl : thickness of warm-layer [m] 107 144 !! 108 145 !! OUTPUT : … … 118 155 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 119 156 !!---------------------------------------------------------------------------------- 157 INTEGER, INTENT(in ) :: kt ! current time step 120 158 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 121 159 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] … … 139 177 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 140 178 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 141 !142 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 180 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 143 181 ! 144 182 INTEGER :: j_itt … … 162 200 !!---------------------------------------------------------------------------------- 163 201 202 IF ( kt == nit000 ) CALL ECMWF_INIT(l_use_cs, l_use_wl) 203 164 204 l_zt_equal_zu = .FALSE. 165 205 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision … … 168 208 IF ( l_use_cs ) THEN 169 209 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN 170 PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!' 171 STOP 172 END IF 173 ALLOCATE ( pdTc(jpi,jpj) ) 174 pdTc(:,:) = -0.25_wp ! First guess of skin correction 210 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP 211 END IF 175 212 END IF 176 213 177 214 IF ( l_use_wl ) THEN 178 IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN 179 PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!' 180 STOP 181 END IF 182 ALLOCATE ( pdTw(jpi,jpj) ) 183 pdTw(:,:) = 0._wp 215 IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN 216 PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP 217 END IF 184 218 END IF 185 219 … … 322 356 !! Cool-skin contribution 323 357 324 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &358 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 325 359 & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0 326 360 327 CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst , pdTc) ! Qnsol -> ztmp1328 329 T_s(:,:) = zsst(:,:) + pdTc(:,:)330 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)361 CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1 362 363 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 364 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 331 365 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 332 366 … … 335 369 IF( l_use_wl ) THEN 336 370 !! Warm-layer contribution 337 CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &371 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 338 372 & ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2 339 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst , pdTw)373 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 340 374 !! Updating T_s and q_s !!! 341 T_s(:,:) = zsst(:,:) + pdTw(:,:)342 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)375 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) 376 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 343 377 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 344 378 END IF 345 346 379 347 380 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN … … 361 394 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 362 395 363 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc 364 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw 396 IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 397 IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 398 IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 365 399 366 400 IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 367 IF ( l_use_cs ) DEALLOCATE ( pdTc )368 IF ( l_use_wl ) DEALLOCATE ( pdTw )369 401 370 402 END SUBROUTINE turb_ecmwf -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
r11638 r11772 45 45 END INTERFACE cp_air 46 46 47 47 INTERFACE alpha_sw 48 48 MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr 49 49 END INTERFACE alpha_sw 50 51 INTERFACE turb_fluxes 52 MODULE PROCEDURE turb_fluxes_vctr, turb_fluxes_sclr 53 END INTERFACE turb_fluxes 50 54 51 55 … … 64 68 PUBLIC update_qnsol_tau 65 69 PUBLIC alpha_sw 70 PUBLIC turb_fluxes 66 71 67 72 !!---------------------------------------------------------------------- … … 128 133 rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 129 134 END FUNCTION rho_air_sclr 130 135 131 136 132 137 … … 464 469 465 470 466 SUBROUTINE UPDATE_QNSOL_TAU( p Ts, pqs, pTa, pqa, pust, ptst, pqst, pUb, pslp, prlw, &471 SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, pslp, prlw, & 467 472 & pQns, pTau, & 468 473 & Qlat) … … 472 477 !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 473 478 !!---------------------------------------------------------------------------------- 479 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 474 480 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 475 481 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 476 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! air temperature at z=zu [K]477 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z= zu [kg/kg]482 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 483 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 478 484 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pust ! u* 479 485 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptst ! t* 480 486 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqst ! q* 481 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=zu [m/s] 487 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 488 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 482 489 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea-level atmospheric pressure [Pa] 483 490 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prlw ! downwelling longwave radiative flux [W/m^2] … … 491 498 & zQlat, zQsen, zQlw 492 499 INTEGER :: ji, jj ! dummy loop indices 493 !!---------------------------------------------------------------------------------- 500 !!---------------------------------------------------------------------------------- 494 501 DO jj = 1, jpj 495 502 DO ji = 1, jpi … … 502 509 zCe = zz0*pqst(ji,jj)/zdq 503 510 504 zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp) ! rho*U10511 !zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp) ! rho*U10 505 512 zTs2 = pTs(ji,jj)*pTs(ji,jj) 506 513 514 CALL TURB_FLUXES( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 515 & pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 516 & pTau(ji,jj), zQsen, zQlat ) 517 518 507 519 ! Wind stress module: 508 pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo?520 !pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 509 521 510 522 ! Non-Solar heat flux to the ocean: 511 zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp ) ! we do not want a Qlat > 0 !512 zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt523 !zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp ) ! we do not want a Qlat > 0 ! 524 !zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt 513 525 zQlw = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 514 526 515 527 pQns(ji,jj) = zQlat + zQsen + zQlw 516 517 IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 528 529 IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 518 530 END DO 519 531 END DO 520 532 END SUBROUTINE UPDATE_QNSOL_TAU 533 534 535 536 537 538 SUBROUTINE TURB_FLUXES_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 539 & pTau, pQsen, pQlat, pEvap ) 540 !!---------------------------------------------------------------------------------- 541 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 542 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 543 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 544 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 545 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 546 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd 547 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCh 548 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCe 549 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 550 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 551 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea-level atmospheric pressure [Pa] 552 !! 553 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2] 554 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen ! [W/m^2] 555 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat ! [W/m^2] 556 !! 557 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! [kg/m^2/s] 558 !! 559 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 560 INTEGER :: ji, jj, jq ! dummy loop indices 561 !!---------------------------------------------------------------------------------- 562 DO jj = 1, jpj 563 DO ji = 1, jpi 564 565 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 566 ztaa = pTa(ji,jj) ! first guess... 567 DO jq = 1, 4 568 zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 569 ztaa = pTa(ji,jj) - zgamma*pzu ! Absolute temp. is slightly colder... 570 END DO 571 zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 572 zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 573 574 zUrho = pUb(ji,jj)*MAX(zrho, 1._wp) ! rho*U10 575 576 pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 577 578 zevap = MIN( zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) , 0._wp ) ! we do not want condensation & Qlat > 0 ! 579 pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 580 pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 581 582 IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 583 584 END DO 585 END DO 586 END SUBROUTINE TURB_FLUXES_VCTR 587 588 589 SUBROUTINE TURB_FLUXES_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 590 & pTau, pQsen, pQlat, pEvap ) 591 !!---------------------------------------------------------------------------------- 592 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 593 REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 594 REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 595 REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 596 REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 597 REAL(wp), INTENT(in) :: pCd 598 REAL(wp), INTENT(in) :: pCh 599 REAL(wp), INTENT(in) :: pCe 600 REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 601 REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 602 REAL(wp), INTENT(in) :: pslp ! sea-level atmospheric pressure [Pa] 603 !! 604 REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2] 605 REAL(wp), INTENT(out) :: pQsen ! [W/m^2] 606 REAL(wp), INTENT(out) :: pQlat ! [W/m^2] 607 !! 608 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! [kg/m^2/s] 609 !! 610 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 611 INTEGER :: jq 612 !!---------------------------------------------------------------------------------- 613 614 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 615 ztaa = pTa ! first guess... 616 DO jq = 1, 4 617 zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) 618 ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder... 619 END DO 620 zrho = rho_air(ztaa, pqa, pslp) 621 zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 622 623 zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10 624 625 pTau = zUrho * pCd * pwnd ! Wind stress module 626 627 zevap = MIN( zUrho * pCe * (pqa - pqs) , 0._wp ) ! we do not want condensation & Qlat > 0 ! 628 pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 629 pQlat = L_vap(pTs) * zevap 630 631 IF ( PRESENT(pEvap) ) pEvap = - zevap 632 633 END SUBROUTINE TURB_FLUXES_SCLR 634 635 636 521 637 522 638 … … 529 645 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 530 646 !!---------------------------------------------------------------------------------- 531 REAL(wp), DIMENSION(jpi,jpj) :: alpha_sw_vctr ! latent heat of vaporization [J/kg]647 REAL(wp), DIMENSION(jpi,jpj) :: alpha_sw_vctr ! thermal expansion coefficient of sea-water [1/K] 532 648 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 533 649 !!---------------------------------------------------------------------------------- … … 543 659 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 544 660 !!---------------------------------------------------------------------------------- 545 REAL(wp) :: alpha_sw_sclr ! latent heat of vaporization [J/kg]661 REAL(wp) :: alpha_sw_sclr ! thermal expansion coefficient of sea-water [1/K] 546 662 REAL(wp), INTENT(in) :: psst ! sea-water temperature [K] 547 663 !!---------------------------------------------------------------------------------- -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
r11672 r11772 38 38 39 39 !! Cool-skin related parameters: 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m] 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 41 & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 42 REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! 43 !! => see eq.(14) in Fairall et al. 1996 (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) 41 44 42 45 !! Warm-layer related parameters: 43 REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp !: maximum depth of warm layer (adjustable) 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 47 & dT_wl, & !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 48 & Hz_wl, & !: depth of warm-layer [m] 49 & Qnt_ac, & !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) 50 & Tau_ac !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) 51 52 REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp !: maximum depth of warm layer (adjustable) 44 53 ! 45 54 REAL(wp), PARAMETER :: rich = 0.65_wp !: critical Richardson number 46 55 ! 47 REAL(wp), PARAMETER :: Qabs_thr = 50._wp !: threshold for heat flux absorbed in WL48 56 REAL(wp), PARAMETER :: zfr0 = 0.5_wp !: initial value of solar flux absorption 49 57 ! 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl ! depth of warm-layer [m]53 58 !!---------------------------------------------------------------------- 54 59 CONTAINS 55 60 56 61 57 SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat , pdT)62 SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat ) 58 63 !!--------------------------------------------------------------------- 59 64 !! 60 65 !! Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019) 61 !! ------------------------------------------------------------------ 66 !! 67 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 68 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 69 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 70 !! doi:10.1029/95JC03190. 71 !! 72 !!------------------------------------------------------------------ 62 73 !! 63 74 !! ** INPUT: … … 68 79 !! *pQlat* surface latent heat flux [K] 69 80 !! 70 !! ** INPUT/OUTPUT:71 !! *pdT* : as input => previous estimate of dT cool-skin72 !! as output => new estimate of dT cool-skin73 !!74 81 !!------------------------------------------------------------------ 75 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 76 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 77 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 78 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 79 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 80 !! 81 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to cool-skin effect 82 !!--------------------------------------------------------------------- 83 INTEGER :: ji, jj ! dummy loop indices 84 REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & 85 & zz1, zz2, zus, & 86 & ztf 87 !!--------------------------------------------------------------------- 88 82 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 83 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 84 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 85 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 86 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 87 !!--------------------------------------------------------------------- 88 INTEGER :: ji, jj, jc 89 REAL(wp) :: zQabs, zdelta, zfr 90 !!--------------------------------------------------------------------- 89 91 DO jj = 1, jpj 90 92 DO ji = 1, jpi 91 93 92 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 93 94 zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere 95 96 zdelta = delta_vl(ji,jj) ! using last value of delta 97 98 !! Fraction of the shortwave flux absorbed by the cool-skin sublayer: 99 zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Eq.16 (Fairall al. 1996b) / !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 100 101 zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 102 103 ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1 104 ! Qt < 0 => warming of the layer => ztf = 0 105 106 !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap): 107 zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0 108 !! LB: this terms only makes sense if > 0 i.e. in the cooling case 109 !! so similar to what's done in ECMWF: 110 zz1 = MAX(0._wp , zz1) ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj)) 111 112 zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases: 113 zz2 = zus*zus * roadrw 114 zz2 = zz2*zz2 115 zlamb = 6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders) 116 117 ! Updating molecular sublayer thickness (delta): 118 zz2 = rnu0_w/(SQRT(roadrw)*zus) 119 zdelta = ztf * zlamb*zz2 & ! Eq.12 (when alpha*Qb>0 / cooling of layer) 120 & + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 ) ! Eq.12 (when alpha*Qb<0 / warming of layer) 121 !LB: changed 0.01 to 0.007 122 123 !! Once again with the new zdelta: 124 zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq.16 (Fairall al. 1996b) 125 zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 126 127 !! Update! 128 pdT(ji,jj) = MIN( - zQnet*zdelta/rk0_w , 0._wp ) ! temperature increment 129 delta_vl(ji,jj) = zdelta 94 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 95 ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 96 97 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 98 99 DO jc = 1, 4 ! because implicit in terms of zdelta... 100 zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) / !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 101 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 102 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 103 END DO 104 105 dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment 130 106 131 107 END DO … … 135 111 136 112 137 138 139 SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait, pdT, Hwl, mask_wl ) 113 SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait ) 140 114 !!--------------------------------------------------------------------- 141 115 !! … … 149 123 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 150 124 !! *iwait* if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop... 151 !! 152 !! ** OUTPUT: 153 !! *pdT* dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST 154 !!--------------------------------------------------------------------- 155 !! 156 !! ** OPTIONAL OUTPUT: 157 !! *Hwl* depth of warm layer [m] 158 !! *mask_wl* mask for possible existence of a warm-layer (1) or not (0) 159 !! 160 !!------------------------------------------------------------------ 125 !!--------------------------------------------------------------------- 161 126 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! surface net solar radiation into the ocean [W/m^2] => >= 0 ! 162 127 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! … … 164 129 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST at depth gdept_1d(1) [K] 165 130 INTEGER , INTENT(in) :: iwait ! if /= 0 then wait before updating accumulated fluxes 166 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST 167 !! 168 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl ! depth of warm layer [m] 169 INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0) 170 ! 171 ! 131 !! 172 132 INTEGER :: ji,jj 173 133 ! 174 REAL(wp) :: zdT _wl, zQabs, zfr, zdz134 REAL(wp) :: zdTwl, zHwl, zQabs, zfr 175 135 REAL(wp) :: zqac, ztac 176 REAL(wp) :: zalpha _w, zcd1, zcd2, flg136 REAL(wp) :: zalpha, zcd1, zcd2, flg 177 137 !!--------------------------------------------------------------------- 178 138 … … 180 140 INTEGER :: jl 181 141 142 LOGICAL :: l_exit, l_destroy_wl 143 182 144 !! INITIALIZATION: 183 pdT(:,:) = 0._wp ! dT initially set to 0._wp 184 zdT_wl = 0._wp ! total warming (amplitude) in warm layer 185 zQabs = 0._wp ! total heat absorped in warm layer 186 zfr = zfr0 ! initial value of solar flux absorption 187 ztac = 0._wp 188 zqac = 0._wp 189 IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0 145 zQabs = 0._wp ! total heat flux absorped in warm layer 146 zfr = zfr0 ! initial value of solar flux absorption !LOLO: save it and use previous value !!! 190 147 191 148 IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize! … … 201 158 DO ji = 1, jpi 202 159 203 zdz = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer 160 l_exit = .FALSE. 161 l_destroy_wl = .FALSE. 162 163 zdTwl = dT_wl(ji,jj) ! value of previous time step as first guess 164 zHwl = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) ! " " " 165 166 zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac 167 ztac = Tau_ac(ji,jj) 204 168 205 169 !***** variables for warm layer *** 206 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 207 208 zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w)) !mess-o-constants 1 209 zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 210 211 !******************************************************** 212 !**** Compute apply warm layer correction ************* 213 !******************************************************** 170 zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 171 172 zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w)) !mess-o-constants 1 173 zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 174 214 175 215 176 znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp ) ! 0<rnoon<1. => rnoon*24 = UTC time of local noon 216 177 zmidn = MOD( znoon-0.5_wp , 1._wp ) 217 218 IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) 'LOLO: rdawn, rdusk, znoon, zmidn =', rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), znoon, zmidn 219 220 ! When midnight is past and dawn is not there yet, do what they call the "midnight reset": 221 IF ( (ztime >= zmidn).AND.(ztime <= rdawn_dcy(ji,jj)) ) THEN 222 IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] MIDNIGHT RESET !!!!, ztime =', ztime 223 zdz = H_wl_max 224 Tau_ac(ji,jj) = 0._wp 225 Qnt_ac(ji,jj) = 0._wp 226 END IF 227 228 229 IF ( ztime > rdawn_dcy(ji,jj) ) THEN ! Dawn, a new day starts at current location 230 IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] WE DO WL !!!!, ztime =', ztime 231 232 !************************************ 233 !**** set warm layer constants *** 234 !************************************ 235 236 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! tot heat absorbed in warm layer 237 238 !PRINT *, ' [WL_COARE] rdt, pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4) 239 240 IF ( zQabs >= Qabs_thr ) THEN ! Check for threshold 241 242 !PRINT *, ' [WL_COARE] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4) 243 244 !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! momentum integral 245 ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral 246 247 IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN !check threshold for warm layer existence 248 !****************************************** 249 ! Compute the absorption profile 250 !****************************************** 251 DO jl = 1, 5 !loop 5 times for zfr 252 zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) & 253 & + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz 254 zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed 255 IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 256 END DO 257 258 ELSE 259 !*********************** 260 ! Warm layer wiped out 261 !*********************** 262 zfr = 0.75 263 zdz = H_wl_max 264 zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed 265 266 END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) 267 268 IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp) !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 ! 269 178 zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 179 180 IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 181 ! Dawn reset to 0! 182 l_exit = .TRUE. 183 l_destroy_wl = .TRUE. 184 END IF 185 186 IF ( .NOT. l_exit ) THEN 187 !! Initial test on initial guess of absorbed heat flux in warm-layer: 188 zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 189 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 190 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!! 191 PRINT *, '#LBD: Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj)) 192 PRINT *, '#LBD: =>Qabs:', zQabs,' zfr=', zfr 193 194 IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 195 ! We have not started to build a WL yet (dT==0) and there's no way it can occur now 196 ! since zQabs <= 0._wp 197 ! => no need to go further 198 PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)' 199 PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs 200 PRINT *, '#LBD: => leaving without changing anything...' 201 l_exit = .TRUE. 202 END IF 203 204 END IF 205 206 ! Okay test on updated absorbed flux: 207 !LOLO: remove??? has a strong influence !!! 208 IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 209 PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt 210 PRINT *, '#LBD: => time to destroy the warm-layer!' 211 l_exit = .TRUE. 212 l_destroy_wl = .TRUE. 213 END IF 214 215 216 IF ( .NOT. l_exit) THEN 217 218 ! Two possibilities at this point: 219 ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0 220 ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it ! 221 222 PRINT *, '#LBD:======================================================' 223 PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4) 224 PRINT *, '#LBD:======================================================' 225 PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4) 226 227 ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral 228 PRINT *, '#LBD: updated value for Tac=', REAL(ztac,4) 229 230 !! We update the value of absorbtion and zQabs: 231 !! some part is useless if Qsw=0 !!! 232 DO jl = 1, 5 233 zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 234 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 235 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) 236 zqac = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed 237 IF ( zqac <= 0._wp ) EXIT 238 zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 239 END DO 240 PRINT *, '#LBD: updated absorption and WL depth=', REAL(zfr,4), REAL(zHwl,4) 241 PRINT *, '#LBD: updated value for Qabs=', REAL(zQabs,4), 'W/m2' 242 PRINT *, '#LBD: updated value for Qac =', REAL(zqac,4), 'J' 243 244 IF ( zqac <= 0._wp ) THEN 245 l_destroy_wl = .TRUE. 246 l_exit = .TRUE. 247 ELSE 248 zdTwl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp) !! => IF(zqac>0._wp): zdTwl=zcd2*zqac**1.5/ztac ; ELSE: zdTwl=0. / ! normally: zqac > 0 ! 249 PRINT *, '#LBD: updated preliminary value for dT_wl=', REAL(zdTwl,4) 270 250 ! Warm layer correction 271 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz) 272 pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 273 274 END IF ! IF ( zQabs >= Qabs_thr ) 275 276 END IF ! IF ( isd_sol >= 21600 ) THEN ! (21600 == 6am) 251 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl) 252 zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) 253 END IF 254 255 END IF !IF ( .NOT. l_exit) 256 257 IF ( l_destroy_wl ) THEN 258 zdTwl = 0._wp 259 zfr = 0.75_wp 260 zHwl = Hwl_max 261 zqac = 0._wp 262 ztac = 0._wp 263 END IF 264 265 PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4) 277 266 278 267 IF ( iwait == 0 ) THEN 279 IF ( (zQabs >= Qabs_thr).AND.(ztime > rdawn_dcy(ji,jj)) ) THEN 280 !PRINT *, ' [WL_COARE] WE UPDATE ACCUMULATED FLUXES !!!' 281 Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 282 Tau_ac(ji,jj) = ztac ! 283 IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1 284 END IF 285 END IF 286 287 H_wl(ji,jj) = zdz 288 289 IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj) 268 !! Iteration loop within bulk algo is over, time to update what needs to be updated: 269 dT_wl(ji,jj) = zdTwl 270 Hz_wl(ji,jj) = zHwl 271 PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4) 272 Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 273 Tau_ac(ji,jj) = ztac 274 PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4) 275 PRINT *, '#LBD' 276 END IF 290 277 291 278 END DO 292 279 END DO 293 280 294 281 END SUBROUTINE WL_COARE 295 282 283 284 285 286 FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a ) 287 !!--------------------------------------------------------------------- 288 !! Computes the thickness (m) of the viscous skin layer. 289 !! Based on Fairall et al., 1996 290 !! 291 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 292 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 293 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 294 !! doi:10.1029/95JC03190. 295 !! 296 !! L. Brodeau, october 2019 297 !!--------------------------------------------------------------------- 298 REAL(wp) :: delta_skin_layer 299 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 300 REAL(wp), INTENT(in) :: pQabs ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 301 REAL(wp), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 302 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 303 !!--------------------------------------------------------------------- 304 REAL(wp) :: zusw, zusw2, zlamb, zQb 305 !!--------------------------------------------------------------------- 306 307 zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay! 308 309 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 310 zusw2 = zusw*zusw 311 312 zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 313 314 delta_skin_layer = zlamb*rnu0_w/zusw 315 316 END FUNCTION delta_skin_layer 296 317 297 318 !!====================================================================== -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
r11637 r11772 38 38 39 39 !! Cool-skin related parameters: 40 ! ... 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 41 & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 42 REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! 43 !! => see eq.(14) in Fairall et al. 1996 (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) 41 44 42 45 !! Warm-layer related parameters: 43 REAL(wp), PARAMETER :: rd0 = 3. !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 44 REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w 45 REAL(wp), PARAMETER :: rNu0 = 1.0 !: be closer to COARE3p6 ???!LOLO 46 !REAL(wp), PARAMETER :: rNu0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 47 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 48 ! !: 0.3 to respect a warming of +3 K in calm 49 ! !: condition for the insolation peak of +1000W/m^2 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 47 & dT_wl, & !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 48 & Hz_wl !: depth of warm-layer [m] 49 ! 50 REAL(wp), PARAMETER, PUBLIC :: rd0 = 3. !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 51 REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w 52 ! 53 REAL(wp), PARAMETER :: rNuwl0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 54 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 55 ! !: 0.3 to respect a warming of +3 K in calm 56 ! !: condition for the insolation peak of +1000W/m^2 50 57 !!---------------------------------------------------------------------- 51 58 CONTAINS 52 59 53 60 54 SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST, pdT ) 55 !!--------------------------------------------------------------------- 56 !! 57 !! Cool-Skin scheme according to Fairall et al. 1996 58 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 59 !! " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 60 !! as parameterized in IFS Cy45r1 / E.C.M.W.F. 61 !! ------------------------------------------------------------------ 61 SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST ) 62 !!--------------------------------------------------------------------- 63 !! 64 !! Cool-skin parameterization, based on Fairall et al., 1996: 65 !! 66 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 67 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 68 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 69 !! doi:10.1029/95JC03190. 70 !! 71 !!------------------------------------------------------------------ 62 72 !! 63 73 !! ** INPUT: … … 66 76 !! *pustar* friction velocity u* [m/s] 67 77 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 68 !!69 !! ** INPUT/OUTPUT:70 !! *pdT* : as input => previous estimate of dT cool-skin71 !! as output => new estimate of dT cool-skin72 !!73 78 !!------------------------------------------------------------------ 74 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 75 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 76 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 77 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 78 !! 79 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to cool-skin effect 80 !!--------------------------------------------------------------------- 81 INTEGER :: ji, jj ! dummy loop indices 82 REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & 83 & zusw, zusw2 84 !!--------------------------------------------------------------------- 85 79 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 80 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 81 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 82 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 83 !!--------------------------------------------------------------------- 84 INTEGER :: ji, jj, jc 85 REAL(wp) :: zQabs, zdelta, zfr 86 !!--------------------------------------------------------------------- 86 87 DO jj = 1, jpj 87 88 DO ji = 1, jpi 88 89 89 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 90 91 zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere 92 93 zusw = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw) ! u* in the water 94 zusw2 = zusw*zusw 95 96 zlamb = 6._wp*( 1._wp + (zQnsol*zalpha_w*rcst_cs/(zusw2*zusw2 ))**0.75 )**(-1./3.) ! w.r.t COARE 3p6 => seems to ommit absorbed zfr*Qsw (Qnet i.o. Qnsol) and effect of evap... 97 ! ! so zlamb not implicit in terms of zdelta (zfr(delta)), so no need to have last guess of delta as in COARE 3.6... 98 zdelta = zlamb*rnu0_w/zusw 99 100 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq. 8.131 / IFS cy40r1, doc, Part IV, 101 zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface 102 103 !! Update! 104 pdT(ji,jj) = MIN( - zQnet*zdelta/rk0_w , 0._wp ) ! temperature increment 90 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 91 ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 92 !zQabs = pQnsol(ji,jj) 93 94 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 95 96 DO jc = 1, 4 ! because implicit in terms of zdelta... 97 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 98 ! => (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 99 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 100 !zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 101 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 102 END DO 103 104 dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment 105 105 106 106 END DO … … 111 111 112 112 113 SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pdT)113 SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pustk ) 114 114 !!--------------------------------------------------------------------- 115 115 !! 116 116 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 117 117 !! " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 118 !! 119 !! STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER! 118 120 !! 119 121 !! As included in IFS Cy45r1 / E.C.M.W.F. … … 125 127 !! *pustar* friction velocity u* [m/s] 126 128 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 127 !!128 !! ** INPUT/OUTPUT:129 !! *pdT* : as input => previous estimate of dT warm-layer130 !! as output => new estimate of dT warm-layer131 !!132 129 !!------------------------------------------------------------------ 133 130 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! surface net solar radiation into the ocean [W/m^2] => >= 0 ! … … 135 132 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity [m/s] 136 133 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST at depth gdept_1d(1) [K] 137 ! 138 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST139 ! 140 INTEGER :: ji, jj134 !! 135 REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pustk ! surface Stokes velocity [m/s] 136 ! 137 INTEGER :: ji, jj, jc 141 138 ! 142 139 REAL(wp) :: & 143 & zdz, & !: thickness of the warm-layer [m] 144 & zalpha_w, & !: thermal expansion coefficient of sea-water 145 & ZSRD, & 146 & dT_wl, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 147 & zfr,zdL,zdL2, ztmp, & 148 & ZPHI, & 149 & zus_a, zusw, zusw2, & 150 & flg, zQabs, ZL1, ZL2 151 !!--------------------------------------------------------------------- 140 & zHwl, & !: thickness of the warm-layer [m] 141 & ztcorr, & !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 142 & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 143 & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 144 & zfr, zeta, ztmp, & 145 & zusw, zusw2, & 146 & zLa, zfLa, & 147 & flg, zwf, zQabs, & 148 & ZA, ZB, zL1, zL2, & 149 & zcst0, zcst1, zcst2, zcst3 150 ! 151 LOGICAL :: l_pustk_known 152 !!--------------------------------------------------------------------- 153 154 l_pustk_known = .FALSE. 155 IF ( PRESENT(pustk) ) l_pustk_known = .TRUE. 152 156 153 157 DO jj = 1, jpj 154 158 DO ji = 1, jpi 155 159 156 zdz = rd0 ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 157 158 ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 160 zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 161 ! it is = rd0 (3m) in default Zeng & Beljaars case... 162 163 !! Previous value of dT / warm-layer, adapted to depth: 164 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 165 ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl 166 zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp ) 167 ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl) 159 168 ! pdT " " and depth of bulk SST (here gdept_1d(1))! 160 !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 161 !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof dT_wl ! 162 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 163 dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 164 165 zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 169 !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced! 170 !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl ! 171 172 zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 166 173 167 174 168 175 ! *** zfr = Fraction of solar radiation absorbed in warm layer (-) 169 zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*z dz) - 0.27_wp*EXP(-2.8_wp*zdz) - 0.45_wp*EXP(-0.07_wp*zdz) !: Eq. 8.157176 zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zHwl) - 0.27_wp*EXP(-2.8_wp*zHwl) - 0.45_wp*EXP(-0.07_wp*zHwl) !: Eq. 8.157 170 177 171 178 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! tot heat absorbed in warm layer 172 179 173 zusw = MAX( pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)! u* in the water180 zusw = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw ! u* in the water 174 181 zusw2 = zusw*zusw 175 182 176 177 !! *** 1st rhs term in eq. 8.156 (IFS doc Cy45r1): 178 ZL1 = zQabs / ( zdz * zRhoCp_w * rNu0 ) * (rNu0 + 1._wp) 179 180 181 !! Buoyancy flux and stability parameter (zdl = -z/L) in water 182 ZSRD = zQabs/zRhoCp_w 183 ! 184 flg = 0.5_wp + SIGN(0.5_wp, ZSRD) ! ZSRD > 0. => 1. / ZSRD < 0. => 0. 185 ztmp = MAX(dT_wl,0._wp) 186 zdl = (flg+1._wp) * ( zusw2 * SQRT(ztmp/(5._wp*zdz*grav*zalpha_w/rNu0)) ) & ! (dT_wl > 0.0 .AND. ZSRD < 0.0) 187 & + flg * ZSRD ! otherwize 188 ! 189 zus_a = MAX( pustar(ji,jj), 1.E-4_wp ) 190 zdL = zdz*vkarmn*grav/(roadrw)**1.5_wp*zalpha_w*zdL/(zus_a*zus_a*zus_a) 191 192 !! Stability function Phi_t(-z/L) (zdL is -z/L) : 193 flg = 0.5_wp + SIGN(0.5_wp, zdL) ! zdl > 0. => 1. / zdl < 0. => 0. 194 zdL2 = zdL*zdL 195 ZPHI = flg * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) & ! (zdL > 0) Takaya et al. 196 & + (flg+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) ) ! (zdl < 0) Eq. 8.136 197 !! FOR zdL > 0.0, old relations: 198 ! ZPHI = 1.+5._wp*zdL ! Eq. 8.136 (Large et al. 1994) 199 ! ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007 200 201 !! *** 2nd rhs term in eq. 8.156 (IFS doc Cy45r1): 202 ZL2 = - (rNu0 + 1._wp) * vkarmn * zusw / ( zdz * ZPHI ) 203 204 ! Forward time / explicit solving of eq. 8.156 (IFS doc Cy45r1): (f_n+1 == pdT(ji,jj) ; f_n == dT_wl) 205 dT_wl = MAX ( dT_wl + rdt*ZL1 + rdt*ZL2*dT_wl , 0._wp ) 206 207 ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz) 208 !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced! 209 210 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$ 211 pdT(ji,jj) = dT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) 212 183 ! Langmuir: 184 IF ( l_pustk_known ) THEN 185 zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6)) 186 ELSE 187 zla = 0.3_wp 188 END IF 189 zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp ) ! Eq.(6) 190 191 zwf = 0.5_wp + SIGN(0.5_wp, zQabs) ! zQabs > 0. => 1. / zQabs < 0. => 0. 192 193 zcst1 = vkarmn*grav*zalpha 194 195 ! 1/L when zQabs > 0 : 196 zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw) 197 198 zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 ) !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ??? 199 200 zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl 201 202 ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w ) 203 204 zcst3 = -zcst0 * vkarmn * zusw * zfLa 205 206 !! Sorry about all these constants ( constant w.r.t zdTwl), it's for 207 !! the sake of optimizations... So all these operations are not done 208 !! over and over within the iteration loop... 209 210 !! T R U L L Y I M P L I C I T => needs itteration 211 !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0.. 212 !! (without this term otherwize the implicit analytical solution is straightforward...) 213 zdTwl_n = zdTwl_b 214 DO jc = 1, 10 215 216 zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence 217 218 ! 1/L when zdTwl > 0 .AND. zQabs < 0 : 219 zL1 = SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw 220 !zL1 = vkarmn*SQRT( zdTwl_n *grav*zalpha / ( 5._wp*zHwl ) ) / zusw ! => vkarmn outside, not inside zcst1 (just for this particular line) ??? 221 222 ! Stability parameter (z/L): 223 zeta = (1._wp - zwf) * zHwl*zL1 + zwf * zHwl*zL2 224 225 ZB = zcst3 / PHI(zeta) 226 227 zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp ) ! Eq.(6) 228 229 END DO 230 231 !! Update: 232 dT_wl(ji,jj) = zdTwl_n * ztcorr 233 213 234 END DO 214 235 END DO … … 216 237 END SUBROUTINE WL_ECMWF 217 238 239 240 241 FUNCTION delta_skin_layer( palpha, pQabs, pustar_a ) 242 !!--------------------------------------------------------------------- 243 !! Computes the thickness (m) of the viscous skin layer. 244 !! Based on Fairall et al., 1996 245 !! 246 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 247 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 248 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 249 !! doi:10.1029/95JC03190. 250 !! 251 !! L. Brodeau, october 2019 252 !!--------------------------------------------------------------------- 253 REAL(wp) :: delta_skin_layer 254 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 255 REAL(wp), INTENT(in) :: pQabs ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 256 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 257 !!--------------------------------------------------------------------- 258 REAL(wp) :: zusw, zusw2, zlamb, zQb 259 !!--------------------------------------------------------------------- 260 261 zQb = pQabs 262 263 !zQ = MIN( -0.1_wp , pQabs ) 264 265 !ztf = 0.5_wp + SIGN(0.5_wp, zQ) ! Qabs < 0 => cooling of the layer => ztf = 0 (normal case) 266 ! ! Qabs > 0 => warming of the layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 267 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 268 zusw2 = zusw*zusw 269 270 zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 271 !zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQ, 0._wp)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 272 273 delta_skin_layer = zlamb*rnu0_w/zusw 274 275 !delta_skin_layer = (1._wp - ztf) * zlamb*rnu0_w/zusw & ! see eq.(12) in Fairall et al., 1996 276 ! & + ztf * MIN(6._wp*rnu0_w/zusw , 0.007_wp) 277 END FUNCTION delta_skin_layer 278 279 280 281 FUNCTION PHI( pzeta) 282 !!--------------------------------------------------------------------- 283 !! 284 !! Takaya et al., 2010 285 !! Eq.(5) 286 !! L. Brodeau, october 2019 287 !!--------------------------------------------------------------------- 288 REAL(wp) :: PHI 289 REAL(wp), INTENT(in) :: pzeta ! stability parameter 290 !!--------------------------------------------------------------------- 291 REAL(wp) :: ztf, zzt2 292 !!--------------------------------------------------------------------- 293 ! 294 zzt2 = pzeta*pzeta 295 ! 296 ztf = 0.5_wp + SIGN(0.5_wp, pzeta) ! zeta > 0 => ztf = 1 297 ! ! zeta < 0 => ztf = 0 298 PHI = ztf * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) & ! zeta > 0 299 & + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) ) ! zeta < 0 300 ! 301 END FUNCTION PHI 302 218 303 !!====================================================================== 219 304 END MODULE sbcblk_skin_ecmwf
Note: See TracChangeset
for help on using the changeset viewer.