Changeset 12015 for NEMO/branches/2019/dev_ASINTER-01-05_merged
- Timestamp:
- 2019-11-29T16:59:07+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_ASINTER-01-05_merged/src
- Files:
-
- 5 added
- 2 deleted
- 8 edited
-
ABL/ablmod.F90 (modified) (5 diffs)
-
ABL/sbcabl.F90 (modified) (1 diff)
-
OCE/SBC/sbc_oce.F90 (modified) (3 diffs)
-
OCE/SBC/sbcblk.F90 (modified) (55 diffs)
-
OCE/SBC/sbcblk_algo_coare.F90 (deleted)
-
OCE/SBC/sbcblk_algo_coare3p0.F90 (added)
-
OCE/SBC/sbcblk_algo_coare3p5.F90 (deleted)
-
OCE/SBC/sbcblk_algo_coare3p6.F90 (added)
-
OCE/SBC/sbcblk_algo_ecmwf.F90 (modified) (13 diffs)
-
OCE/SBC/sbcblk_algo_ncar.F90 (modified) (17 diffs)
-
OCE/SBC/sbcblk_phy.F90 (added)
-
OCE/SBC/sbcblk_skin_coare.F90 (added)
-
OCE/SBC/sbcblk_skin_ecmwf.F90 (added)
-
OCE/SBC/sbcdcy.F90 (modified) (8 diffs)
-
OCE/SBC/sbcmod.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ABL/ablmod.F90
r11937 r12015 17 17 USE phycst ! physical constants 18 18 USE dom_oce, ONLY : tmask 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 20 USE sbcblk ! use some physical constants for flux computation 19 USE sbc_oce, ONLY : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 20 USE sbcblk ! use rn_?fac 21 USE sbcblk_phy ! use some physical constants for flux computation 21 22 ! 22 23 USE prtctl ! Print control (prt_ctl routine) … … 100 101 REAL(wp) , INTENT( out), DIMENSION(:,: ) :: ptauj_ice ! ice-surface tauy stress (V-point) 101 102 #endif 102 103 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: z rhoa, zwnd_i, zwnd_j103 ! 104 REAL(wp), DIMENSION(1:jpi,1:jpj ) :: zwnd_i, zwnd_j 104 105 REAL(wp), DIMENSION(1:jpi,2:jpka ) :: zCF 105 106 REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka) :: z_cft !--FL--to be removed after the test phase … … 529 530 ztemp = tq_abl ( ji, jj, 2, nt_a, jp_ta ) 530 531 zhumi = tq_abl ( ji, jj, 2, nt_a, jp_qa ) 531 zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 532 & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 532 !zcff = pslp_dta( ji, jj ) / & !<-- At this point ztemp and zhumi should not be zero ... 533 ! & ( R_dry*ztemp * ( 1._wp + rctv0*zhumi ) ) 534 zcff = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 533 535 psen ( ji, jj ) = cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 534 536 pevp ( ji, jj ) = rn_efac*MAX( 0._wp, zcff * pevp(ji,jj) * ( pssq(ji,jj) - zhumi ) ) 535 zrhoa( ji, jj ) = zcff537 rhoa( ji, jj ) = zcff 536 538 END DO 537 539 END DO … … 551 553 zcff = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) & 552 554 & + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) ! * msk_abl(ji,jj) 553 zztmp = zrhoa(ji,jj) * pcd_du(ji,jj)555 zztmp = rhoa(ji,jj) * pcd_du(ji,jj) 554 556 555 557 pwndm (ji,jj) = zcff … … 593 595 zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 594 596 595 ptaui_ice(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) &596 & + zrhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) &597 ptaui_ice(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj) & 598 & + rhoa(ji ,jj) * pCd_du_ice(ji ,jj) ) & 597 599 & * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 598 ptauj_ice(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) &599 & + zrhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) &600 ptauj_ice(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1) & 601 & + rhoa(ji,jj ) * pCd_du_ice(ji,jj ) ) & 600 602 & * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 601 603 END DO -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ABL/sbcabl.F90
r11858 r12015 341 341 & tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), & ! <<= in 342 342 & sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m , & ! <<= in 343 & sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) , & ! <<= in 343 344 & zssq, zcd_du, zsen, zevp ) ! =>> out 344 345 -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbc_oce.F90
r11275 r12015 116 116 !! wndm is used compute surface gases exchanges in ice-free ocean or leads 117 117 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wndm !: wind speed module at T-point (=|U10m-Uoce|) [m/s] 118 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhoa !: air density at "rn_zu" m above the sea [kg/m3] !LB 118 119 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qsr !: sea heat flux: solar [W/m2] 119 120 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qns , qns_b !: sea heat flux: non solar [W/m2] … … 159 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_m !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 160 161 162 !!---------------------------------------------------------------------- 163 !! Cool-skin/Warm-layer 164 !!---------------------------------------------------------------------- 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tsk !: sea-surface skin temperature (used if ln_skin_cs==.true. .OR. ln_skin_wl==.true.) [K] !LB 166 167 161 168 !! * Substitutions 162 169 # include "vectopt_loop_substitute.h90" … … 177 184 ! 178 185 ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) , & 179 & vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) )186 & vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) ) 180 187 ! 181 188 ALLOCATE( qns_tot(jpi,jpj) , qns (jpi,jpj) , qns_b(jpi,jpj), & -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90
r11586 r12015 15 15 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 16 16 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 17 !! ! ==> based on AeroBulk (http ://aerobulk.sourceforge.net/)17 !! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 20 20 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 21 21 !!---------------------------------------------------------------------- … … 24 24 !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition 25 25 !! sbc_blk : bulk formulation as ocean surface boundary condition 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 28 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 29 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 30 !! q_sat : saturation humidity as a function of SLP and temperature 31 !! L_vap : latent heat of vaporization of water as a function of temperature 32 !! sea-ice case only : 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 28 !! sea-ice case only : 33 29 !! blk_ice_1 : provide the air-ice stress 34 30 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 35 31 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 36 32 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 37 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 33 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 38 34 !!---------------------------------------------------------------------- 39 35 USE oce ! ocean dynamics and tracers … … 51 47 USE icethd_dh ! for CALL ice_thd_snwblow 52 48 #endif 53 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 54 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)55 USE sbcblk_algo_coare3p 5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)56 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 50 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 51 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 52 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1) 57 53 ! 58 54 USE iom ! I/O manager library … … 62 58 USE prtctl ! Print control 63 59 60 USE sbcblk_phy !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc... 61 62 64 63 IMPLICIT NONE 65 64 PRIVATE … … 69 68 PUBLIC blk_oce_1 ! called in sbcabl 70 69 PUBLIC blk_oce_2 ! called in sbcabl 71 PUBLIC rho_air ! called in ablmod72 PUBLIC cp_air ! called in ablmod73 70 #if defined key_si3 74 71 PUBLIC blk_ice_1 ! routine called in icesbc 75 72 PUBLIC blk_ice_2 ! routine called in icesbc 76 73 PUBLIC blk_ice_qcn ! routine called in icesbc 77 #endif 78 79 INTERFACE cp_air 80 MODULE PROCEDURE cp_air_0d, cp_air_2d 81 END INTERFACE 82 83 ! !!* Namelist namsbc_blk : bulk parameters 84 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 85 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 86 LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013) 87 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31) 88 ! ! 89 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 90 REAL(wp), PUBLIC :: rn_efac !: multiplication factor for evaporation 91 REAL(wp), PUBLIC :: rn_vfac !: multiplication factor for ice/ocean velocity in the calculation of wind stress 92 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 93 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 94 ! ! 95 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 96 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 97 98 INTEGER :: nblk ! choice of the bulk algorithm 99 ! ! associated indices: 100 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 101 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 102 INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013) 103 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31) 104 105 ! !!! air parameters 106 REAL(wp) , PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg] 107 REAL(wp) , PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg] 108 REAL(wp), PUBLIC, PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 109 REAL(wp) , PARAMETER :: R_vap = 461.495_wp ! Specific gas constant for water vapor [J/K/kg] 110 REAL(wp) , PARAMETER :: reps0 = R_dry/R_vap ! ratio of gas constant for dry air and water vapor => ~ 0.622 111 REAL(wp), PUBLIC, PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 112 ! !!! Bulk parameters 113 REAL(wp) , PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...) 114 REAL(wp) , PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 115 REAL(wp) , PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 116 REAL(wp) , PARAMETER :: rcdice = 1.4e-3 ! transfer coefficient over ice 117 REAL(wp) , PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 74 #endif 75 76 INTEGER , PUBLIC :: jpfld ! maximum number of files to read 77 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 79 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 80 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 81 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 83 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 85 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 86 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 88 89 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 90 91 ! !!* Namelist namsbc_blk : bulk parameters 92 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 93 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 94 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 95 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 96 ! 97 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 98 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 99 REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 100 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 101 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 102 ! 103 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 104 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 118 105 ! 119 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice … … 121 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 122 109 123 !INTEGER , PUBLIC, PARAMETER :: jpfld =11 !: maximum number of files to read 124 INTEGER , PUBLIC :: jpfld !: maximum number of files to read 125 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 !: index of 10m wind velocity (i-component) (m/s) at T-point 126 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 !: index of 10m wind velocity (j-component) (m/s) at T-point 127 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 !: index of 10m air temperature (Kelvin) 128 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 !: index of specific humidity ( % ) 129 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 !: index of solar heat (W/m2) 130 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 !: index of Long wave (W/m2) 131 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 !: index of total precipitation (rain+snow) (Kg/m2/s) 132 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 !: index of snow (solid prcipitation) (kg/m2/s) 133 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 !: index of sea level pressure (Pa) 134 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 135 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 136 137 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 110 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 111 LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 112 LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 113 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 114 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 115 ! 116 INTEGER :: nhumi ! choice of the bulk algorithm 117 ! ! associated indices: 118 INTEGER, PARAMETER :: np_humi_sph = 1 119 INTEGER, PARAMETER :: np_humi_dpt = 2 120 INTEGER, PARAMETER :: np_humi_rlh = 3 121 122 INTEGER :: nblk ! choice of the bulk algorithm 123 ! ! associated indices: 124 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 125 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 126 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 127 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1) 138 128 139 129 !! * Substitutions … … 165 155 !! ** Purpose : choose and initialize a bulk formulae formulation 166 156 !! 167 !! ** Method : 157 !! ** Method : 168 158 !! 169 159 !!---------------------------------------------------------------------- … … 172 162 !! 173 163 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 174 !TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 175 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 164 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 176 165 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 177 166 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " … … 179 168 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 180 169 & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & 181 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 182 & cn_dir , rn_zqt, rn_zu, & 183 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 170 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 171 & cn_dir , rn_zqt, rn_zu, & 172 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 173 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 184 174 !!--------------------------------------------------------------------- 185 175 ! … … 187 177 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 188 178 ! 189 ! !** read bulk namelist 179 ! !** read bulk namelist 190 180 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 191 181 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) … … 200 190 ! !** initialization of the chosen bulk formulae (+ check) 201 191 ! !* select the bulk chosen in the namelist and check the choice 202 ioptio = 0 203 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF 204 IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF 205 IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF 206 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 207 ! 192 ioptio = 0 193 IF( ln_NCAR ) THEN 194 nblk = np_NCAR ; ioptio = ioptio + 1 195 ENDIF 196 IF( ln_COARE_3p0 ) THEN 197 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 198 ENDIF 199 IF( ln_COARE_3p6 ) THEN 200 nblk = np_COARE_3p6 ; ioptio = ioptio + 1 201 ENDIF 202 IF( ln_ECMWF ) THEN 203 nblk = np_ECMWF ; ioptio = ioptio + 1 204 ENDIF 208 205 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 206 207 ! !** initialization of the cool-skin / warm-layer parametrization 208 IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) & 209 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 210 ! 211 ioptio = 0 212 IF( ln_humi_sph ) THEN 213 nhumi = np_humi_sph ; ioptio = ioptio + 1 214 ENDIF 215 IF( ln_humi_dpt ) THEN 216 nhumi = np_humi_dpt ; ioptio = ioptio + 1 217 ENDIF 218 IF( ln_humi_rlh ) THEN 219 nhumi = np_humi_rlh ; ioptio = ioptio + 1 220 ENDIF 221 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 209 222 ! 210 223 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 211 224 IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 212 IF( sn_qsr%ln_tint ) THEN 225 IF( sn_qsr%ln_tint ) THEN 213 226 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 214 227 & ' ==> We force time interpolation = .false. for qsr' ) … … 234 247 ALLOCATE( sf(jpfld), STAT=ierror ) 235 248 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 236 ! !- fill the bulk structure with namelist informations237 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )238 249 ! 239 250 DO jfpr= 1, jpfld … … 258 269 ENDIF 259 270 END DO 260 ! 261 IF( ln_wave ) THEN ! surface waves 262 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) & ! Activated wave module but neither drag nor stokes drift activated 263 & CALL ctl_stop( 'sbc_blk_init: Ask for wave coupling but ln_cdgw=ln_sdw=ln_tauwoc=ln_stcor=F' ) 264 IF( ln_cdgw .AND. .NOT.ln_NCAR ) & ! drag coefficient read from wave model only with NCAR bulk formulae 265 & CALL ctl_stop( 'sbc_blk_init: drag coefficient read from wave model need NCAR bulk formulae') 266 IF( ln_stcor .AND. .NOT.ln_sdw ) & 267 CALL ctl_stop( 'sbc_blk_init: Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 271 ! !- fill the bulk structure with namelist informations 272 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 273 ! 274 IF( ln_wave ) THEN 275 !Activated wave module but neither drag nor stokes drift activated 276 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 277 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 278 !drag coefficient read from wave model definable only with mfs bulk formulae and core 279 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 280 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 281 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 282 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 283 ENDIF 268 284 ELSE 269 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) THEN 270 CALL ctl_warn( 'sbc_blk_init: ln_wave=F, set all wave-related namelist parameter to FALSE') 271 ln_cdgw =.FALSE. ; ln_sdw =.FALSE. ; ln_tauwoc =.FALSE. ; ln_stcor =.FALSE. 272 ENDIF 273 ENDIF 285 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 286 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 287 & 'with drag coefficient (ln_cdgw =T) ' , & 288 & 'or Stokes Drift (ln_sdw=T) ' , & 289 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 290 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 291 ENDIF 274 292 ! 275 293 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient … … 281 299 ! 282 300 ! set transfer coefficients to default sea-ice values 283 Cd_ice(:,:) = r cdice284 Ch_ice(:,:) = r cdice285 Ce_ice(:,:) = r cdice301 Cd_ice(:,:) = rCd_ice 302 Ch_ice(:,:) = rCd_ice 303 Ce_ice(:,:) = rCd_ice 286 304 ! 287 305 IF(lwp) THEN !** Control print 288 306 ! 289 WRITE(numout,*) !* namelist 307 WRITE(numout,*) !* namelist 290 308 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 291 309 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 292 310 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 293 WRITE(numout,*) ' "COARE 3. 5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0294 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31)ln_ECMWF = ', ln_ECMWF311 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 312 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 295 313 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 296 314 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 306 324 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 307 325 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 308 CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)' 309 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)' 326 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 327 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 328 END SELECT 329 ! 330 WRITE(numout,*) 331 WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs 332 WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl 333 ! 334 WRITE(numout,*) 335 SELECT CASE( nhumi ) !* Print the choice of air humidity 336 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 337 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 338 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 310 339 END SELECT 311 340 ! … … 355 384 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 356 385 ! 357 ! ! compute the surface ocean fluxes using bulk formulea358 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN359 !386 IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1) ! no previous estimate of skin temperature => using bulk SST 387 ! 388 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 360 389 CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in 361 390 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in 362 391 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 392 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 363 393 & zssq, zcd_du, zsen, zevp ) ! =>> out 364 394 … … 368 398 & zsen, zevp ) ! <=> in out 369 399 ENDIF 370 400 ! 371 401 #if defined key_cice 372 IF( MOD( kt -1, nn_fsbc ) == 0 ) THEN402 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 373 403 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 374 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 375 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 376 ENDIF 404 IF( ln_dm2dc ) THEN 405 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 406 ELSE 407 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 408 ENDIF 377 409 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 378 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 410 411 SELECT CASE( nhumi ) 412 CASE( np_humi_sph ) 413 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 414 CASE( np_humi_dpt ) 415 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 416 CASE( np_humi_rlh ) 417 qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 418 END SELECT 419 379 420 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 380 421 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 387 428 388 429 389 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp 390 & pslp, pst , pu , pv, & ! inp 391 & pssq, pcd_du, psen, pevp ) ! out 430 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp 431 & pslp , pst , pu , pv, & ! inp 432 & pqsr , pqlw , & ! inp 433 & pssq , pcd_du, psen , pevp ) ! out 392 434 !!--------------------------------------------------------------------- 393 435 !! *** ROUTINE blk_oce_1 *** 394 436 !! 395 !! ** Purpose : Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 396 !! 397 !! ** Method : bulk formulae using atmospheric 398 !! fields from the ABL model at previous time-step 437 !! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes 438 !! if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 439 !! 440 !! ** Method : bulk formulae using atmospheric fields from : 441 !! if ln_blk=T, atmospheric fields read in sbc_read 442 !! if ln_abl=T, the ABL model at previous time-step 399 443 !! 400 444 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) … … 412 456 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 413 457 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 458 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 459 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 414 460 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 415 461 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] … … 423 469 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 424 470 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 425 REAL(wp), DIMENSION(jpi,jpj) :: z rhoa ! density of air [kg/m^3]471 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] 426 472 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 427 473 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 428 474 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 475 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux 476 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 429 477 !!--------------------------------------------------------------------- 430 478 ! … … 460 508 461 509 ! ----------------------------------------------------------------------------- ! 462 ! I Turbulent FLUXES!510 ! I Solar FLUX ! 463 511 ! ----------------------------------------------------------------------------- ! 464 512 465 ! ... specific humidity at SST and IST tmask( 466 pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) ) 467 ! 513 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 514 zztmp = 1. - albo 515 IF( ln_dm2dc ) THEN 516 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 517 ELSE 518 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 519 ENDIF 520 521 522 ! ----------------------------------------------------------------------------- ! 523 ! II Turbulent FLUXES ! 524 ! ----------------------------------------------------------------------------- ! 525 526 ! specific humidity at SST 527 pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 528 529 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 530 zztmp1(:,:) = zst(:,:) 531 zztmp2(:,:) = pssq(:,:) 532 ENDIF 533 534 ! specific humidity of air at "rn_zqt" m above the sea 535 SELECT CASE( nhumi ) 536 CASE( np_humi_sph ) 537 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 538 CASE( np_humi_dpt ) 539 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 540 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 541 CASE( np_humi_rlh ) 542 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 543 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 544 END SELECT 545 ! 546 ! potential temperature of air at "rn_zqt" m above the sea 468 547 IF( ln_abl ) THEN 469 548 ztpot = ptair(:,:) … … 472 551 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 473 552 ! (since reanalysis products provide T at z, not theta !) 474 ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt 475 ENDIF 476 477 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 478 ! 479 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! NCAR-COREv2 480 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 481 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.0 482 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 483 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.5 484 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 485 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! ECMWF 486 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 487 CASE DEFAULT 488 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 489 END SELECT 490 491 492 IF ( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 493 !! FL do we need this multiplication by tmask ... ??? 494 DO jj = 1, jpj 495 DO ji = 1, jpi 496 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 497 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 498 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 499 psen(ji,jj) = zztmp * zch_oce(ji,jj) 500 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 501 END DO 502 END DO 503 504 ELSE !== BLK formulation ==! turbulent fluxes computation 553 ztpot = ptair(:,:) + gamma_moist( ptair(:,:), zqair(:,:) ) * rn_zqt 554 ENDIF 555 556 557 558 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 559 560 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 561 562 CASE( np_COARE_3p0 ) 563 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 564 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 565 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 566 567 CASE( np_COARE_3p6 ) 568 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 569 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 570 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 571 572 CASE( np_ECMWF ) 573 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 574 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 575 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 576 577 CASE DEFAULT 578 CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' ) 579 END SELECT 580 581 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 582 WHERE ( fr_i < 0.001_wp ) 583 ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 584 zst(:,:) = zst(:,:)*tmask(:,:,1) 585 pssq(:,:) = pssq(:,:)*tmask(:,:,1) 586 ELSEWHERE 587 ! we forget about the update... 588 zst(:,:) = zztmp1(:,:) !LB: using what we backed up before skin-algo 589 pssq(:,:) = zztmp2(:,:) !LB: " " " 590 END WHERE 591 592 !LB: Update of tsk, the "official" array for skin temperature 593 tsk(:,:) = zst(:,:) 594 595 ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 596 597 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 598 ! 599 CASE( np_NCAR ) 600 CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, & 601 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 602 603 CASE( np_COARE_3p0 ) 604 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 605 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 606 607 CASE( np_COARE_3p6 ) 608 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 609 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 610 611 CASE( np_ECMWF ) 612 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 613 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 614 615 CASE DEFAULT 616 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 617 END SELECT 618 619 ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 620 621 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef. 622 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef. 623 624 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 625 !! If zu == zt, then ensuring once for all that: 626 t_zu(:,:) = ztpot(:,:) 627 q_zu(:,:) = zqair(:,:) 628 ENDIF 629 630 631 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 632 ! ------------------------------------------------------------- 633 634 IF (ln_abl) THEN 635 636 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 637 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 638 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 639 & pcd_du(:,:), psen(:,:), pevp(:,:), & 640 & prhoa=rhoa(:,:) ) 641 642 ELSE 643 644 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 645 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 646 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 647 & taum(:,:), psen(:,:), zqla(:,:), & 648 & pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 649 650 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 651 psen(:,:) = psen(:,:) * tmask(:,:,1) 652 taum(:,:) = taum(:,:) * tmask(:,:,1) 653 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 505 654 506 ! ! Compute true air density : 507 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 508 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) ) 509 ELSE ! At zt: 510 zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 511 END IF 512 513 DO jj = 1, jpj 514 DO ji = 1, jpi 515 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) 516 !!gm to be done by someone: check the efficiency of the call of cp_air within do loops 517 psen (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * zch_oce(ji,jj) * ( zst(ji,jj) - t_zu(ji,jj) ) 518 pevp (ji,jj) = rn_efac*MAX( 0._wp, zztmp * zce_oce(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) ) 519 zztmp = zztmp * zcd_oce(ji,jj) 520 pcd_du(ji,jj) = zztmp 521 taum (ji,jj) = zztmp * wndm (ji,jj) 522 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 523 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 524 END DO 525 END DO 655 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 656 zcd_oce = 0._wp 657 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 658 zwnd_i = zcd_oce * zwnd_i 659 zwnd_j = zcd_oce * zwnd_j 526 660 527 661 CALL iom_put( "taum_oce", taum ) ! output wind stress module 528 662 529 663 ! ... utau, vtau at U- and V_points, resp. 530 664 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines … … 539 673 END DO 540 674 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 541 542 675 IF(ln_ctl) THEN 543 676 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') … … 554 687 ENDIF 555 688 ! 556 END SUBROUTINE blk_oce_1689 END SUBROUTINE blk_oce_1 557 690 558 691 … … 593 726 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 594 727 728 595 729 ! ----------------------------------------------------------------------------- ! 596 ! I Radiative FLUXES!730 ! III Net longwave radiative FLUX ! 597 731 ! ----------------------------------------------------------------------------- ! 598 732 599 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 600 zztmp = 1._wp - albo 601 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1) 602 ELSE ; qsr(:,:) = zztmp * pqsr(:,:) * tmask(:,:,1) 603 ENDIF 604 605 zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 733 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 734 !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 735 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 606 736 607 737 ! Turbulent fluxes over ocean 608 738 ! ----------------------------- 609 739 610 zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux 740 zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 611 741 612 742 IF(ln_ctl) THEN 613 743 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 614 744 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 745 615 746 ENDIF 616 747 617 748 ! ----------------------------------------------------------------------------- ! 618 ! I IITotal FLUXES !749 ! IV Total FLUXES ! 619 750 ! ----------------------------------------------------------------------------- ! 620 751 ! 621 emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1) ! mass flux (evap. - precip.) 622 ! 623 zz1 = rn_pfac * rLfus 624 zz2 = rn_pfac * rcp 625 zz3 = rn_pfac * rcpi 626 zztmp = 1._wp - albo 627 qns(:,:) = ( zqlw(:,:) - psen(:,:) - zqla(:,:) & ! Downward Non Solar 628 & - psnow(:,:) * zztmp & ! remove latent melting heat for solid precip 629 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 630 & + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2 & ! add liquid precip heat content at Tair 631 & + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3 & ! add solid precip heat content at min(Tair,Tsnow) 632 & ) * tmask(:,:,1) 752 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 753 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 754 ! 755 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 756 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 757 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? 758 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 759 & * ( ptair(:,:) - rt0 ) * rcp & 760 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 761 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 762 qns(:,:) = qns(:,:) * tmask(:,:,1) 633 763 ! 634 764 #if defined key_si3 635 qns_oce(:,:) = zqlw(:,:) - psen(:,:) -zqla(:,:) ! non solar without emp (only needed by SI3)765 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 636 766 qsr_oce(:,:) = qsr(:,:) 637 767 #endif 638 768 ! 639 IF ( nn_ice == 0 ) THEN 640 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 641 CALL iom_put( "qsb_oce" , - psen ) ! output downward sensible heat over the ocean 642 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 643 CALL iom_put( "qemp_oce", qns-zqlw+psen+zqla ) ! output downward heat content of E-P over the ocean 644 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 645 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 646 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 647 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 648 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 649 CALL iom_put( 'snowpre', sprecip ) ! Snow 650 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 769 CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB 770 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 771 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 772 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 773 CALL iom_put( "evap_oce" , pevp ) ! evaporation 774 CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean 775 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 776 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 777 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 778 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 779 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 780 CALL iom_put( 'snowpre', sprecip ) ! Snow 781 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 782 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 783 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 784 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 651 785 ENDIF 652 786 ! … … 660 794 661 795 662 FUNCTION rho_air( ptak, pqa, pslp )663 !!-------------------------------------------------------------------------------664 !! *** FUNCTION rho_air ***665 !!666 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere667 !!668 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)669 !!-------------------------------------------------------------------------------670 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]671 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]672 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]673 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]674 !!-------------------------------------------------------------------------------675 !676 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) )677 !678 END FUNCTION rho_air679 680 681 FUNCTION cp_air_0d( pqa )682 !!-------------------------------------------------------------------------------683 !! *** FUNCTION cp_air ***684 !!685 !! ** Purpose : Compute specific heat (Cp) of moist air686 !!687 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)688 !!-------------------------------------------------------------------------------689 REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]690 REAL(wp) :: cp_air_0d! specific heat of moist air [J/K/kg]691 !!-------------------------------------------------------------------------------692 !693 Cp_air_0d = Cp_dry + Cp_vap * pqa694 !695 END FUNCTION cp_air_0d696 697 698 FUNCTION cp_air_2d( pqa )699 !!-------------------------------------------------------------------------------700 !! *** FUNCTION cp_air ***701 !!702 !! ** Purpose : Compute specific heat (Cp) of moist air703 !!704 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)705 !!-------------------------------------------------------------------------------706 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]707 REAL(wp), DIMENSION(jpi,jpj) :: cp_air_2d! specific heat of moist air [J/K/kg]708 !!-------------------------------------------------------------------------------709 !710 Cp_air_2d = Cp_dry + Cp_vap * pqa711 !712 END FUNCTION cp_air_2d713 714 715 FUNCTION q_sat( ptak, pslp )716 !!----------------------------------------------------------------------------------717 !! *** FUNCTION q_sat ***718 !!719 !! ** Purpose : Specific humidity at saturation in [kg/kg]720 !! Based on accurate estimate of "e_sat"721 !! aka saturation water vapor (Goff, 1957)722 !!723 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)724 !!----------------------------------------------------------------------------------725 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]726 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]727 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]728 !729 INTEGER :: ji, jj ! dummy loop indices730 REAL(wp) :: ze_sat, ztmp ! local scalar731 !!----------------------------------------------------------------------------------732 !733 DO jj = 1, jpj734 DO ji = 1, jpi735 !736 ztmp = rt0 / ptak(ji,jj)737 !738 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)739 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &740 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) &741 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 )742 !743 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa]744 !745 END DO746 END DO747 !748 END FUNCTION q_sat749 750 751 FUNCTION gamma_moist( ptak, pqa )752 !!----------------------------------------------------------------------------------753 !! *** FUNCTION gamma_moist ***754 !!755 !! ** Purpose : Compute the moist adiabatic lapse-rate.756 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate757 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html758 !!759 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)760 !!----------------------------------------------------------------------------------761 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]762 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]763 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate764 !765 INTEGER :: ji, jj ! dummy loop indices766 REAL(wp) :: zrv, ziRT ! local scalar767 !!----------------------------------------------------------------------------------768 !769 DO jj = 1, jpj770 DO ji = 1, jpi771 zrv = pqa(ji,jj) / (1. - pqa(ji,jj))772 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT773 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )774 END DO775 END DO776 !777 END FUNCTION gamma_moist778 779 780 FUNCTION L_vap( psst )781 !!---------------------------------------------------------------------------------782 !! *** FUNCTION L_vap ***783 !!784 !! ** Purpose : Compute the latent heat of vaporization of water from temperature785 !!786 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)787 !!----------------------------------------------------------------------------------788 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]789 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]790 !!----------------------------------------------------------------------------------791 !792 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6793 !794 END FUNCTION L_vap795 796 796 #if defined key_si3 797 797 !!---------------------------------------------------------------------- … … 802 802 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 803 803 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 804 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 804 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 805 805 !!---------------------------------------------------------------------- 806 806 … … 835 835 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 836 836 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 837 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau)838 837 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 839 838 !!--------------------------------------------------------------------- 839 ! 840 840 841 841 ! ------------------------------------------------------------ ! … … 858 858 Ce_ice(:,:) = Cd_ice(:,:) 859 859 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 860 CALL Cdn10_Lupkes2015( ptsui, Cd_ice, Ch_ice )860 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 861 861 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 862 862 ENDIF 863 863 864 IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef.865 IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef.864 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 865 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 866 866 867 867 ! local scalars ( place there for vector optimisation purposes) 868 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 869 zrhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 868 IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 870 869 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 871 870 … … 877 876 DO jj = 2, jpjm1 878 877 DO ji = fs_2, fs_jpim1 ! vect. opt. 879 putaui(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * zcd_dui(ji+1,jj) &880 & + zrhoa(ji ,jj) * zcd_dui(ji ,jj) ) &878 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 879 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 881 880 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 882 pvtaui(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * zcd_dui(ji,jj+1) &883 & + zrhoa(ji,jj ) * zcd_dui(ji,jj ) ) &881 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 882 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 884 883 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 885 884 END DO … … 898 897 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 899 898 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 900 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj)899 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 901 900 END DO 902 901 END DO … … 934 933 REAL(wp) :: zst3 ! local variable 935 934 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 936 REAL(wp) :: zztmp, zztmp2, z1_rLsub 935 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 937 936 REAL(wp) :: zfr1, zfr2 ! local variables 938 937 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 942 941 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 943 942 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 944 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa 945 !!--------------------------------------------------------------------- 946 ! 947 zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars 948 zcoef_dqla = -Ls * 11637800. * (-5897.8) 949 ! 950 zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 943 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 944 !!--------------------------------------------------------------------- 945 ! 946 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 947 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 948 ! 949 SELECT CASE( nhumi ) 950 CASE( np_humi_sph ) 951 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 952 CASE( np_humi_dpt ) 953 zqair(:,:) = q_sat( phumi(:,:), pslp ) 954 CASE( np_humi_rlh ) 955 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 956 END SELECT 951 957 ! 952 958 zztmp = 1. / ( 1. - albo ) 953 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 954 ELSEWHERE ; z1_st(:,:,:) = 0._wp 959 WHERE( ptsu(:,:,:) /= 0._wp ) 960 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 961 ELSEWHERE 962 z1_st(:,:,:) = 0._wp 955 963 END WHERE 956 964 ! ! ========================== ! … … 966 974 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 967 975 ! Long Wave (lw) 968 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)976 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 969 977 ! lw sensitivity 970 978 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 976 984 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 977 985 ! Sensible Heat 978 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))986 z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 979 987 ! Latent Heat 980 988 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 981 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls* Ce_ice(ji,jj) * wndm_ice(ji,jj) * &982 & ( 11637800. * zztmp2 / zrhoa(ji,jj) - phumi(ji,jj) ) )989 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 990 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 983 991 ! Latent heat sensitivity for ice (Dqla/Dt) 984 992 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN … … 990 998 991 999 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 992 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_ice(ji,jj) * wndm_ice(ji,jj)993 1000 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 1001 994 1002 ! ----------------------------! 995 1003 ! III Total FLUXES ! … … 1042 1050 DO jl = 1, jpl 1043 1051 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 1044 ! ! But we do not have Tice => consider it at 0degC => evap=0 1052 ! ! But we do not have Tice => consider it at 0degC => evap=0 1045 1053 END DO 1046 1054 … … 1049 1057 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 1050 1058 ! 1051 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1059 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1052 1060 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 1053 1061 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 1054 1062 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 1055 1063 ELSEWHERE ! zero when hs>0 1056 qtr_ice_top(:,:,:) = 0._wp 1064 qtr_ice_top(:,:,:) = 0._wp 1057 1065 END WHERE 1058 1066 ! … … 1064 1072 CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl) 1065 1073 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 1066 END 1074 ENDIF 1067 1075 ! 1068 1076 END SUBROUTINE blk_ice_2 1069 1077 1070 1078 1071 1079 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 1076 1084 !! to force sea ice / snow thermodynamics 1077 1085 !! in the case conduction flux is emulated 1078 !! 1086 !! 1079 1087 !! ** Method : compute surface energy balance assuming neglecting heat storage 1080 1088 !! following the 0-layer Semtner (1976) approach … … 1101 1109 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 1102 1110 !!--------------------------------------------------------------------- 1103 1111 1104 1112 ! -------------------------------------! 1105 1113 ! I Enhanced conduction factor ! … … 1109 1117 ! 1110 1118 zgfac(:,:,:) = 1._wp 1111 1119 1112 1120 IF( ld_virtual_itd ) THEN 1113 1121 ! … … 1115 1123 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 1116 1124 zfac3 = 2._wp / zepsilon 1117 ! 1118 DO jl = 1, jpl 1125 ! 1126 DO jl = 1, jpl 1119 1127 DO jj = 1 , jpj 1120 1128 DO ji = 1, jpi … … 1124 1132 END DO 1125 1133 END DO 1126 ! 1127 ENDIF 1128 1134 ! 1135 ENDIF 1136 1129 1137 ! -------------------------------------------------------------! 1130 1138 ! II Surface temperature and conduction flux ! … … 1136 1144 DO jj = 1 , jpj 1137 1145 DO ji = 1, jpi 1138 ! 1146 ! 1139 1147 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1140 1148 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) … … 1153 1161 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1154 1162 qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) & 1155 1163 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1156 1164 1157 1165 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1158 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1166 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1159 1167 1160 1168 END DO 1161 1169 END DO 1162 1170 ! 1163 END DO 1164 ! 1171 END DO 1172 ! 1165 1173 END SUBROUTINE blk_ice_qcn 1166 1174 1167 1175 1168 1176 SUBROUTINE Cdn10_Lupkes2012( pcd ) … … 1170 1178 !! *** ROUTINE Cdn10_Lupkes2012 *** 1171 1179 !! 1172 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1180 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1173 1181 !! to make it dependent on edges at leads, melt ponds and flows. 1174 1182 !! After some approximations, this can be resumed to a dependency 1175 1183 !! on ice concentration. 1176 !! 1184 !! 1177 1185 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1178 1186 !! with the highest level of approximation: level4, eq.(59) … … 1186 1194 !! 1187 1195 !! This new drag has a parabolic shape (as a function of A) starting at 1188 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1196 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1189 1197 !! and going down to Cdi(say 1.4e-3) for A=1 1190 1198 !! … … 1211 1219 1212 1220 ! ice-atm drag 1213 pcd(:,:) = r cdice + & ! pure ice drag1221 pcd(:,:) = rCd_ice + & ! pure ice drag 1214 1222 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1215 1223 1216 1224 END SUBROUTINE Cdn10_Lupkes2012 1217 1225 1218 1226 1219 SUBROUTINE Cdn10_Lupkes2015( ptm_su, p cd, pch )1227 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 1220 1228 !!---------------------------------------------------------------------- 1221 1229 !! *** ROUTINE Cdn10_Lupkes2015 *** 1222 1230 !! 1223 1231 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1224 !! between sea-ice and atmosphere with distinct momentum 1225 !! and heat coefficients depending on sea-ice concentration 1232 !! between sea-ice and atmosphere with distinct momentum 1233 !! and heat coefficients depending on sea-ice concentration 1226 1234 !! and atmospheric stability (no meltponds effect for now). 1227 !! 1235 !! 1228 1236 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1229 1237 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1230 1238 !! it considers specific skin and form drags (Andreas et al. 2010) 1231 !! to compute neutral transfert coefficients for both heat and 1239 !! to compute neutral transfert coefficients for both heat and 1232 1240 !! momemtum fluxes. Atmospheric stability effect on transfert 1233 1241 !! coefficient is also taken into account following Louis (1979). … … 1239 1247 ! 1240 1248 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] 1249 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa] 1241 1250 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient 1242 1251 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient … … 1267 1276 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw 1268 1277 REAL(wp) :: zCdn_form_tmp 1269 !!---------------------------------------------------------------------- -----1270 ! 1278 !!---------------------------------------------------------------------- 1279 1271 1280 ! Momentum Neutral Transfert Coefficients (should be a constant) 1272 1281 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 … … 1277 1286 ! Heat Neutral Transfert Coefficients 1278 1287 zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 1279 1288 1280 1289 ! Atmospheric and Surface Variables 1281 1290 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1282 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ocean [kg/kg]1283 zqi_sat(:,:) = 0.98_wp * q_sat( ptm_su(:,:), sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ice [kg/kg]1284 ! 1285 DO jj = 2, jpjm1 ! reduced loop is necessary for reproduc tibility1291 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg] 1292 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1293 ! 1294 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1286 1295 DO ji = fs_2, fs_jpim1 1287 1296 ! Virtual potential temperature [K] … … 1289 1298 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1290 1299 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1291 1300 1292 1301 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1293 1302 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1294 1303 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1295 1304 1296 1305 ! Momentum and Heat Neutral Transfert Coefficients 1297 1306 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1298 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1299 1300 ! Momentum and Heat Stability functions ( possibility to use psi_m_ecmwf instead)1307 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1308 1309 ! Momentum and Heat Stability functions (!!GS: possibility to use psi_m_ecmwf instead ?) 1301 1310 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1302 1311 z0i = z0_skin_ice ! over ice … … 1309 1318 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1310 1319 ENDIF 1311 1320 1312 1321 IF( zrib_i <= 0._wp ) THEN 1313 1322 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 … … 1317 1326 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1318 1327 ENDIF 1319 1328 1320 1329 ! Momentum Transfert Coefficients (Eq. 38) 1321 1330 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1322 1331 & zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1323 1332 1324 1333 ! Heat Transfert Coefficients (Eq. 49) 1325 1334 pch(ji,jj) = zChn_skin_ice * zfhi + & -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r10069 r12015 1 1 MODULE sbcblk_algo_ecmwf 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes turbulent components of surface fluxes 5 !! according to the method in IFS of the ECMWF model 6 !! 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2)10 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 13 11 !! based on IFS doc (avaible online on the ECMWF's website) 14 12 !! 13 !! Routine turb_ecmwf maintained and developed in AeroBulk 14 !! (https://github.com/brodeau/aerobulk) 15 15 !! 16 !! Routine turb_ecmwf maintained and developed in AeroBulk 17 !! (http://aerobulk.sourceforge.net/) 18 !! 19 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 16 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 20 17 !!---------------------------------------------------------------------- 21 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 41 38 42 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 41 USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 43 42 44 43 IMPLICIT NONE 45 44 PRIVATE 46 45 47 PUBLIC :: TURB_ECMWF ! called by sbcblk.F9048 49 ! !! ECMWF own values for given constants, taken form IFS documentation...46 PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 47 48 !! ECMWF own values for given constants, taken form IFS documentation... 50 49 REAL(wp), PARAMETER :: charn0 = 0.018 ! Charnock constant (pretty high value here !!! 51 50 ! ! => Usually 0.011 for moderate winds) 52 51 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 53 52 REAL(wp), PARAMETER :: Beta0 = 1. ! gustiness parameter ( = 1.25 in COAREv3) 54 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...55 REAL(wp), PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg]56 REAL(wp), PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg]57 53 REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term) 58 54 REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1) 59 55 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 56 57 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 58 60 59 !!---------------------------------------------------------------------- 61 60 CONTAINS 62 61 63 SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu, & 64 & Cd, Ch, Ce , t_zu, q_zu, U_blk, & 65 & Cdn, Chn, Cen ) 66 !!---------------------------------------------------------------------------------- 67 !! *** ROUTINE turb_ecmwf *** 68 !! 69 !! 2015: L. Brodeau (brodeau@gmail.com) 70 !! 71 !! ** Purpose : Computes turbulent transfert coefficients of surface 72 !! fluxes according to IFS doc. (cycle 31) 73 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 74 !! 75 !! ** Method : Monin Obukhov Similarity Theory 62 63 SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 64 !!--------------------------------------------------------------------- 65 !! *** FUNCTION sbcblk_algo_ecmwf_init *** 76 66 !! 77 67 !! INPUT : 78 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 ) CALL ctl_stop( ' SBCBLK_ALGO_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 ENDIF 83 IF( l_use_cs ) THEN 84 ierr = 0 85 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 86 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 87 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 88 ENDIF 89 END SUBROUTINE sbcblk_algo_ecmwf_init 90 91 92 93 SUBROUTINE turb_ecmwf( 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, & 95 & Cdn, Chn, Cen, & 96 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 97 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 98 !!---------------------------------------------------------------------- 99 !! *** ROUTINE turb_ecmwf *** 100 !! 101 !! ** Purpose : Computes turbulent transfert coefficients of surface 102 !! fluxes according to IFS doc. (cycle 45r1) 103 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 104 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 105 !! 106 !! Applies the cool-skin warm-layer correction of the SST to T_s 107 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 108 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 109 !! are provided as (optional) arguments! 110 !! 111 !! INPUT : 112 !! ------- 113 !! * kt : current time step (starts at 1) 79 114 !! * zt : height for temperature and spec. hum. of air [m] 80 !! * zu : height for wind speed (generally 10m) [m] 81 !! * U_zu : scalar wind speed at 10m [m/s] 82 !! * sst : SST [K] 115 !! * zu : height for wind speed (usually 10m) [m] 83 116 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 117 !! * q_zt : specific humidity of air at zt [kg/kg] 86 !! 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 121 !! 122 !! INPUT/OUTPUT: 123 !! ------------- 124 !! * T_s : always "bulk SST" as input [K] 125 !! -> unchanged "bulk SST" as output if CSWL not used [K] 126 !! -> skin temperature as output if CSWL used [K] 127 !! 128 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 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: 133 !! --------------- 134 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 135 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 136 !! * slp : sea-level pressure [Pa] 137 !! 138 !! OPTIONAL OUTPUT: 139 !! ---------------- 140 !! * pdT_cs : SST increment "dT" for cool-skin correction [K] 141 !! * pdT_wl : SST increment "dT" for warm-layer correction [K] 142 !! * pHz_wl : thickness of warm-layer [m] 87 143 !! 88 144 !! OUTPUT : … … 93 149 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 150 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m [m/s] 96 !! 97 !! 98 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 99 !!---------------------------------------------------------------------------------- 151 !! * U_blk : bulk wind speed at zu [m/s] 152 !! 153 !! 154 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 155 !!---------------------------------------------------------------------------------- 156 INTEGER, INTENT(in ) :: kt ! current time step 100 157 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 158 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]159 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 160 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq! sea surface specific humidity [kg/kg]105 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity 161 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 162 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 163 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 164 LOGICAL , INTENT(in ) :: l_use_cs ! use the cool-skin parameterization 165 LOGICAL , INTENT(in ) :: l_use_wl ! use the warm-layer parameterization 107 166 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 108 167 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) … … 110 169 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 170 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 112 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind at 10m[m/s]171 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 172 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 173 ! 174 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 175 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 176 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 177 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_cs 178 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pdT_wl ! [K] 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 180 ! 115 181 INTEGER :: j_itt 116 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 118 ! 119 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star, & 120 & dt_zu, dq_zu, & 121 & znu_a, & !: Nu_air, Viscosity of air 122 & Linv, & !: 1/L (inverse of Monin Obukhov length... 123 & z0, z0t, z0q 124 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 125 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 126 !!---------------------------------------------------------------------------------- 127 ! 128 ! Identical first gess as in COARE, with IFS parameter values though 129 ! 182 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 183 ! 184 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 185 REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 186 REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 187 REAL(wp), DIMENSION(jpi,jpj) :: Linv !: 1/L (inverse of Monin Obukhov length... 188 REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 189 ! 190 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst ! to back up the initial bulk SST 191 ! 192 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 193 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 194 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 195 !!---------------------------------------------------------------------------------- 196 197 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 198 130 199 l_zt_equal_zu = .FALSE. 131 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 132 133 200 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 201 202 !! Initializations for cool skin and warm layer: 203 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 204 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 205 206 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 207 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 208 209 IF( l_use_cs .OR. l_use_wl ) THEN 210 ALLOCATE ( zsst(jpi,jpj) ) 211 zsst = T_s ! backing up the bulk SST 212 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 213 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 214 ENDIF 215 216 217 ! Identical first gess as in COARE, with IFS parameter values though... 218 ! 134 219 !! First guess of temperature and humidity at height zu: 135 t_zu = MAX( t_zt , 0.0) ! who knows what's given on masked-continental regions...136 q_zu = MAX( q_zt , 1.e-6 ) ! "220 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 221 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 137 222 138 223 !! Pot. temp. difference (and we don't want it to be 0!) 139 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu )140 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu )141 142 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)143 144 ztmp2 = 0.5 * 0.5! initial guess for wind gustiness contribution145 U_blk = SQRT(U_zu*U_zu + ztmp2) 146 147 ! z0 = 0.0001148 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0149 ztmp0 = LOG(zu*ztmp2) 150 z tmp1 = LOG(10.*ztmp2)151 u_star = 0.035*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10)152 153 z0 = charn0*u_star*u_star/grav + 0.11*znu_a/u_star154 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t !224 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 225 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 226 227 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 228 229 U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 230 231 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 232 ztmp1 = LOG(10._wp*10000._wp) ! " " " 233 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 234 235 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 238 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 239 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 155 240 156 241 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 157 242 158 ztmp0 = vkarmn*vkarmn/LOG(zt *z0t)/Cd159 160 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number161 162 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu(ztmp2):163 ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )243 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 244 245 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 246 247 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 248 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 164 249 func_m = ztmp0*ztmp2 ! temporary array !! 165 !! Ribu < 0 Ribu > 0 Beta = 1.25166 func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! temporary array !!! func_h == zeta_u167 & + ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0))250 func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u 251 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 252 !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 168 253 169 254 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 170 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h))171 172 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h))255 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 256 257 u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 173 258 t_star = dt_zu*ztmp0 174 259 q_star = dq_zu*ztmp0 175 260 176 ! What 's needto be done if zt /= zu:261 ! What needs to be done if zt /= zu: 177 262 IF( .NOT. l_zt_equal_zu ) THEN 178 !179 263 !! First update of values at zu (or zt for wind) 180 264 ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu) ! zt*func_h/zu == zeta_t 181 ztmp1 = log(zt/zu) + ztmp0265 ztmp1 = LOG(zt/zu) + ztmp0 182 266 t_zu = t_zt - t_star/vkarmn*ztmp1 183 267 q_zu = q_zt - q_star/vkarmn*ztmp1 184 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 185 186 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 187 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 268 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 188 269 ! 270 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 271 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 189 272 ENDIF 190 273 … … 194 277 195 278 !! First guess of inverse of Monin-Obukov length (1/L) : 196 ztmp0 = (1. + rctv0*q_zu) ! the factor to apply to temp. to get virt. temp... 197 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 279 Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 198 280 199 281 !! Functions such as u* = U_blk*vkarmn/func_m 200 ztmp1 = zu + z0 201 ztmp0 = ztmp1*Linv 202 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 203 func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 204 282 ztmp0 = zu*Linv 283 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 284 func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 205 285 206 286 !! ITERATION BLOCK 207 !! ***************208 209 287 DO j_itt = 1, nb_itt 210 288 211 289 !! Bulk Richardson Number at z=zu (Eq. 3.25) 212 ztmp0 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)290 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 213 291 214 292 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 215 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 293 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 294 !! Note: it is slightly different that the L we would get with the usual 295 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 216 296 217 297 !! Update func_m with new Linv: 218 ztmp1 = zu + z0 219 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 298 func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 220 299 221 300 !! Need to update roughness lengthes: … … 223 302 ztmp2 = u_star*u_star 224 303 ztmp1 = znu_a/u_star 225 z0 = alpha_M*ztmp1 + charn0*ztmp2/grav 226 z0t = alpha_H*ztmp1 ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 227 z0q = alpha_Q*ztmp1 228 229 !! Update wind at 10m taking into acount convection-related wind gustiness: 230 ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 231 ztmp2 = ztmp2 * (MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 232 !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 233 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 305 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 306 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 307 308 !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 309 ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 310 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 311 U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 234 312 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 235 313 … … 238 316 !! as well the air-sea differences: 239 317 IF( .NOT. l_zt_equal_zu ) THEN 240 241 318 !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 242 func_h = psi_h_ecmwf( (zu+z0)*Linv) ! temporary array !!!243 func_m = psi_h_ecmwf( (zt+z0)*Linv) ! temporary array !!!319 func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 320 func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 244 321 245 322 ztmp2 = psi_h_ecmwf(z0t*Linv) 246 323 ztmp0 = func_h - ztmp2 247 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0t) - ztmp0)324 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 248 325 t_star = dt_zu*ztmp1 249 326 ztmp2 = ztmp0 - func_m + ztmp2 … … 253 330 ztmp2 = psi_h_ecmwf(z0q*Linv) 254 331 ztmp0 = func_h - ztmp2 255 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0q) - ztmp0)332 ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 256 333 q_star = dq_zu*ztmp1 257 334 ztmp2 = ztmp0 - func_m + ztmp2 258 ztmp1 = log(zt/zu) + ztmp2335 ztmp1 = LOG(zt/zu) + ztmp2 259 336 q_zu = q_zt - q_star/vkarmn*ztmp1 260 261 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 262 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 263 264 END IF 337 ENDIF 265 338 266 339 !! Updating because of updated z0 and z0t and new Linv... 267 ztmp1 = zu + z0 268 ztmp0 = ztmp1*Linv 269 func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 270 func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 271 272 END DO 340 ztmp0 = zu*Linv 341 func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 342 func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 343 344 345 IF( l_use_cs ) THEN 346 !! Cool-skin contribution 347 348 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, & 349 & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0 350 351 CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst ) ! Qnsol -> ztmp1 352 353 T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 354 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 355 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 356 357 ENDIF 358 359 IF( l_use_wl ) THEN 360 !! Warm-layer contribution 361 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, & 362 & ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2 363 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 364 !! Updating T_s and q_s !!! 365 T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) ! 366 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 367 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 368 ENDIF 369 370 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 371 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 372 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 373 ENDIF 374 375 END DO !DO j_itt = 1, nb_itt 273 376 274 377 Cd = vkarmn*vkarmn/(func_m*func_m) 275 378 Ch = vkarmn*vkarmn/(func_m*func_h) 276 ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 277 Ce = vkarmn*vkarmn/(func_m*ztmp1) 278 279 ztmp1 = zu + z0 280 Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 281 Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 282 Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 283 284 END SUBROUTINE TURB_ECMWF 379 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 380 Ce = vkarmn*vkarmn/(func_m*ztmp2) 381 382 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 383 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 384 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 385 386 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 387 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 388 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 389 390 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 391 392 END SUBROUTINE turb_ecmwf 285 393 286 394 … … 294 402 !! and L is M-O length 295 403 !! 296 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)404 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 297 405 !!---------------------------------------------------------------------------------- 298 406 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 302 410 REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 303 411 !!---------------------------------------------------------------------------------- 304 !305 412 DO jj = 1, jpj 306 413 DO ji = 1, jpi 307 414 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):415 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 309 416 ! 310 417 ! Unstable (Paulson 1970): 311 418 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta))313 ztmp = 1. + SQRT(zx)419 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 420 ztmp = 1._wp + SQRT(zx) 314 421 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125 *ztmp*(1.+ zx) ) &316 & -2. *ATAN( SQRT(zx) ) + 0.5*rpi422 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 423 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 317 424 ! 318 425 ! Unstable: 319 426 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 320 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) &321 & - zzeta - 2. /3.*5./0.35427 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 428 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 322 429 ! 323 430 ! Combining: 324 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1325 ! 326 psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable327 & + stab * psi_stab ! (zzeta > 0) Stable431 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 432 ! 433 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 434 & + stab * psi_stab ! (zzeta > 0) Stable 328 435 ! 329 436 END DO 330 437 END DO 331 !332 438 END FUNCTION psi_m_ecmwf 333 439 334 440 335 441 FUNCTION psi_h_ecmwf( pzeta ) 336 442 !!---------------------------------------------------------------------------------- … … 342 448 !! and L is M-O length 343 449 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)450 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 451 !!---------------------------------------------------------------------------------- 346 452 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 354 460 DO ji = 1, jpi 355 461 ! 356 zzeta = MIN(pzeta(ji,jj) , 5. ) ! Very stable conditions (L positif and big!):357 ! 358 zx = ABS(1. - 16.*zzeta)**.25 ! this is actually (1/phi_m)**2 !!!462 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 463 ! 464 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 359 465 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 466 ! Unstable (Paulson 1970) : 361 psi_unst = 2. *LOG(0.5*(1.+ zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1467 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 468 ! 363 469 ! Stable: 364 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1365 & - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.470 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 471 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 366 472 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 473 ! 368 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1369 ! 370 ! 371 psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable372 & + stab * psi_stab ! (zzeta > 0) Stable474 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 475 ! 476 ! 477 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 478 & + stab * psi_stab ! (zzeta > 0) Stable 373 479 ! 374 480 END DO 375 481 END DO 376 !377 482 END FUNCTION psi_h_ecmwf 378 483 379 380 FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )381 !!----------------------------------------------------------------------------------382 !! Bulk Richardson number (Eq. 3.25 IFS doc)383 !!384 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)385 !!----------------------------------------------------------------------------------386 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk !387 !388 REAL(wp) , INTENT(in) :: pz ! height above the sea [m]389 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptz ! air temperature at pz m [K]390 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdt ! ptz - sst [K]391 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqz ! air temperature at pz m [kg/kg]392 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdq ! pqz - ssq [kg/kg]393 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s]394 !!----------------------------------------------------------------------------------395 !396 Ri_bulk = grav*pz/(pub*pub) &397 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz))) &398 & + rctv0*pdq )399 !400 END FUNCTION Ri_bulk401 402 403 FUNCTION visc_air(ptak)404 !!----------------------------------------------------------------------------------405 !! Air kinetic viscosity (m^2/s) given from temperature in degrees...406 !!407 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)408 !!----------------------------------------------------------------------------------409 REAL(wp), DIMENSION(jpi,jpj) :: visc_air !410 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)411 !412 INTEGER :: ji, jj ! dummy loop indices413 REAL(wp) :: ztc, ztc2 ! local scalar414 !!----------------------------------------------------------------------------------415 !416 DO jj = 1, jpj417 DO ji = 1, jpi418 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C419 ztc2 = ztc*ztc420 visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)421 END DO422 END DO423 !424 END FUNCTION visc_air425 484 426 485 !!====================================================================== -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ncar.F90
r10190 r12015 11 11 !! 12 12 !! Routine turb_ncar maintained and developed in AeroBulk 13 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk/) 14 14 !! 15 15 !! L. Brodeau, 2015 … … 26 26 USE dom_oce ! ocean space and time domain 27 27 USE phycst ! physical constants 28 USE sbc_oce ! Surface boundary condition: ocean fields 28 USE iom ! I/O manager library 29 USE lib_mpp ! distribued memory computing library 30 USE in_out_manager ! I/O manager 31 USE prtctl ! Print control 29 32 USE sbcwave, ONLY : cdn_wave ! wave module 30 33 #if defined key_si3 || defined key_cice 31 34 USE sbc_ice ! Surface boundary condition: ice fields 32 35 #endif 33 !34 USE iom ! I/O manager library35 USE lib_mpp ! distribued memory computing library36 USE in_out_manager ! I/O manager37 USE prtctl ! Print control38 36 USE lib_fortran ! to use key_nosignedzero 39 37 38 USE sbc_oce ! Surface boundary condition: ocean fields 39 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 40 41 41 IMPLICIT NONE 42 42 PRIVATE 43 43 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 46 ! ! NCAR own values for given constants: 47 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature... 48 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 46 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 47 49 48 !!---------------------------------------------------------------------- 50 49 CONTAINS … … 53 52 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 54 53 & Cdn, Chn, Cen ) 55 !!---------------------------------------------------------------------- ------------54 !!---------------------------------------------------------------------- 56 55 !! *** ROUTINE turb_ncar *** 57 56 !! … … 61 60 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 62 61 !! 63 !! ** Method : Monin Obukhov Similarity Theory64 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)65 !!66 !! ** References : Large & Yeager, 2004 / Large & Yeager, 200867 !!68 !! ** Last update: Laurent Brodeau, June 2014:69 !! - handles both cases zt=zu and zt/=zu70 !! - optimized: less 2D arrays allocated and less operations71 !! - better first guess of stability by checking air-sea difference of virtual temperature72 !! rather than temperature difference only...73 !! - added function "cd_neutral_10m" that uses the improved parametrization of74 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!75 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them76 !! => 'vkarmn' and 'grav'77 !!78 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)79 62 !! 80 63 !! INPUT : 81 64 !! ------- 82 65 !! * zt : height for temperature and spec. hum. of air [m] 83 !! * zu : height for wind speed (generally 10m) [m] 84 !! * U_zu : scalar wind speed at 10m [m/s] 85 !! * sst : SST [K] 66 !! * zu : height for wind speed (usually 10m) [m] 67 !! * sst : bulk SST [K] 86 68 !! * t_zt : potential air temperature at zt [K] 87 69 !! * ssq : specific humidity at saturation at SST [kg/kg] 88 70 !! * q_zt : specific humidity of air at zt [kg/kg] 71 !! * U_zu : scalar wind speed at zu [m/s] 89 72 !! 90 73 !! … … 96 79 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 97 80 !! * q_zu : specific humidity of air // [kg/kg] 98 !! * U_blk : bulk wind at 10m [m/s] 81 !! * U_blk : bulk wind speed at zu [m/s] 82 !! 83 !! 84 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 85 !!---------------------------------------------------------------------------------- 100 86 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 103 89 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 90 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq ! sea surface specific humidity [kg/kg] 105 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity 91 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 92 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 93 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 110 96 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 111 97 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 112 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind at 10m[m/s]98 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 113 99 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 100 ! 115 INTEGER :: j_itt 116 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 101 INTEGER :: j_itt 102 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 118 103 ! 119 104 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10 ! 10m neutral latent/sensible coefficient … … 126 111 ! 127 112 l_zt_equal_zu = .FALSE. 128 IF( ABS(zu - zt) < 0.01 )l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision129 130 U_blk = MAX( 0.5 , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s113 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 114 115 U_blk = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 131 116 132 117 !! First guess of stability: 133 ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt134 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable118 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 119 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 135 120 136 121 !! Neutral coefficients at 10m: … … 139 124 ztmp0 (:,:) = cdn_wave(:,:) 140 125 ELSE 141 126 ztmp0 = cd_neutral_10m( U_blk ) 142 127 ENDIF 143 128 … … 146 131 !! Initializing transf. coeff. with their first guess neutral equivalents : 147 132 Cd = ztmp0 148 Ce = 1.e-3 *( 34.6* sqrt_Cd_n10 )149 Ch = 1.e-3 *sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab))133 Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 ) 134 Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) 150 135 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd) 151 136 152 IF( ln_cdgw ) Cen = Ce ; Chn = Ch 137 IF( ln_cdgw ) THEN 138 Cen = Ce 139 Chn = Ch 140 ENDIF 153 141 154 142 !! Initializing values at z_u with z_t values: 155 143 t_zu = t_zt ; q_zu = q_zt 156 144 157 !! * Now starting iteration loop158 DO j_itt =1, nb_itt145 !! ITERATION BLOCK 146 DO j_itt = 1, nb_itt 159 147 ! 160 148 ztmp1 = t_zu - sst ! Updating air/sea differences … … 162 150 163 151 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 164 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 165 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 166 167 ztmp0 = 1. + rctv0*q_zu ! multiply this with t and you have the virtual temperature 152 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd)) 153 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 154 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 168 155 169 156 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 170 ztmp0 = (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 171 ! ( Cd*U_blk*U_blk is U*^2 at zu ) 172 157 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 158 173 159 !! Stability parameters : 174 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 160 zeta_u = zu*ztmp0 161 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 175 162 zpsi_h_u = psi_h( zeta_u ) 176 163 … … 178 165 IF( .NOT. l_zt_equal_zu ) THEN 179 166 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 180 stab = zt*ztmp0 ; stab = SIGN( MIN(ABS(stab),10.0), stab ) ! Temporaty array stab == zeta_t !!! 167 stab = zt*ztmp0 168 stab = SIGN( MIN(ABS(stab),10._wp), stab ) ! Temporaty array stab == zeta_t !!! 181 169 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 182 170 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 183 171 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 184 q_zu = max(0., q_zu) 185 END IF 186 172 q_zu = max(0._wp, q_zu) 173 ENDIF 174 175 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 176 ! In very rare low-wind conditions, the old way of estimating the 177 ! neutral wind speed at 10m leads to a negative value that causes the code 178 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 187 179 ztmp2 = psi_m(zeta_u) 188 180 IF( ln_cdgw ) THEN ! surface wave case 189 181 stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd)) 190 182 Cd = stab * stab 191 ztmp0 = (LOG(zu/10. ) - zpsi_h_u) / vkarmn / sqrt_Cd_n10183 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 192 184 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 193 ztmp1 = 1. + Chn * ztmp0185 ztmp1 = 1._wp + Chn * ztmp0 194 186 Ch = Chn * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 195 ztmp1 = 1. + Cen * ztmp0187 ztmp1 = 1._wp + Cen * ztmp0 196 188 Ce = Cen * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 197 189 198 190 ELSE 199 200 201 202 203 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))204 205 206 207 208 stab = 0.5 + sign(0.5,zeta_u)! update stability209 Cx_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.- stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10)210 211 212 213 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))214 215 216 217 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10218 219 ztmp1 = 1.+ Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10)220 221 222 Cx_n10 = 1.e-3 * (34.6* sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10223 224 ztmp1 = 1.+ Cx_n10*ztmp0225 226 227 ! 228 END DO 229 ! 191 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 192 ! In very rare low-wind conditions, the old way of estimating the 193 ! neutral wind speed at 10m leads to a negative value that causes the code 194 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 195 ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 196 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 197 Cdn(:,:) = ztmp0 198 sqrt_Cd_n10 = sqrt(ztmp0) 199 200 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 201 Cx_n10 = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 202 Chn(:,:) = Cx_n10 203 204 !! Update of transfer coefficients: 205 ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 206 Cd = ztmp0 / ( ztmp1*ztmp1 ) 207 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 208 209 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 210 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 211 ztmp1 = 1._wp + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 212 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 213 214 Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 215 Cen(:,:) = Cx_n10 216 ztmp1 = 1._wp + Cx_n10*ztmp0 217 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 218 ENDIF 219 220 END DO !DO j_itt = 1, nb_itt 221 230 222 END SUBROUTINE turb_ncar 231 223 … … 238 230 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 239 231 !! 240 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)232 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 241 233 !!---------------------------------------------------------------------------------- 242 234 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) … … 255 247 ! 256 248 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 257 zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) ) ! If pw10 < 33. => 0, else => 1258 ! 259 cd_neutral_10m(ji,jj) = 1.e-3 * ( &260 & (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind < 33 m/s261 & + zgt33 * 2.34 )! wind >= 33 m/s262 ! 263 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6 )249 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 250 ! 251 cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 252 & (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind < 33 m/s 253 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s 254 ! 255 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 264 256 ! 265 257 END DO … … 273 265 !! Universal profile stability function for momentum 274 266 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 275 !! 276 !! pzet 0 : stability paramenter, z/L where z is altitude measurement267 !! 268 !! pzeta : stability paramenter, z/L where z is altitude measurement 277 269 !! and L is M-O length 278 270 !! 279 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)280 !!---------------------------------------------------------------------------------- 281 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pzeta282 REAL(wp), DIMENSION(jpi,jpj) :: psi_m283 ! 284 INTEGER :: ji, jj 271 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 272 !!---------------------------------------------------------------------------------- 273 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 274 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 275 ! 276 INTEGER :: ji, jj ! dummy loop indices 285 277 REAL(wp) :: zx2, zx, zstab ! local scalars 286 278 !!---------------------------------------------------------------------------------- 287 !288 279 DO jj = 1, jpj 289 280 DO ji = 1, jpi 290 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )291 zx2 = MAX ( zx2 , 1.)281 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 282 zx2 = MAX( zx2 , 1._wp ) 292 283 zx = SQRT( zx2 ) 293 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )294 ! 295 psi_m(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable296 & + (1. - zstab) * (2.*LOG((1. + zx)*0.5) & ! Unstable297 & + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5) ! "284 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 285 ! 286 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 287 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 288 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 298 289 ! 299 290 END DO 300 291 END DO 301 !302 292 END FUNCTION psi_m 303 293 … … 308 298 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 309 299 !! 310 !! pzet 0 : stability paramenter, z/L where z is altitude measurement300 !! pzeta : stability paramenter, z/L where z is altitude measurement 311 301 !! and L is M-O length 312 302 !! 313 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 314 !!---------------------------------------------------------------------------------- 303 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 304 !!---------------------------------------------------------------------------------- 305 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 315 306 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 316 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 307 ! 308 INTEGER :: ji, jj ! dummy loop indices 319 309 REAL(wp) :: zx2, zstab ! local scalars 320 310 !!---------------------------------------------------------------------------------- … … 322 312 DO jj = 1, jpj 323 313 DO ji = 1, jpi 324 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )325 zx2 = MAX ( zx2 , 1.)326 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )327 ! 328 psi_h(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable329 & + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5)) ! Unstable314 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 315 zx2 = MAX( zx2 , 1._wp ) 316 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 317 ! 318 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 319 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 330 320 ! 331 321 END DO 332 322 END DO 333 !334 323 END FUNCTION psi_h 335 324 -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcdcy.F90
r10425 r12015 7 7 !! NEMO 2.0 ! 2006-02 (S. Masson, G. Madec) adaptation to NEMO 8 8 !! 3.1 ! 2009-07 (J.M. Molines) adaptation to v3.1 9 !! 4.* ! 2019-10 (L. Brodeau) nothing really new, but the routine 10 !! ! "sbc_dcy_param" has been extracted from old function "sbc_dcy" 11 !! ! => this allows the warm-layer param of COARE3* to know the time 12 !! ! of dawn and dusk even if "ln_dm2dc=.false." (rdawn_dcy & rdusk_dcy 13 !! ! are now public) 9 14 !!---------------------------------------------------------------------- 10 15 … … 22 27 IMPLICIT NONE 23 28 PRIVATE 24 29 25 30 INTEGER, PUBLIC :: nday_qsr !: day when parameters were computed 26 31 27 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: raa , rbb , rcc , rab ! diurnal cycle parameters 28 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtmd, rdawn, rdusk, rscal ! - - - 29 33 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtmd, rscal ! - - - 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rdawn_dcy, rdusk_dcy ! - - - 35 30 36 PUBLIC sbc_dcy ! routine called by sbc 37 PUBLIC sbc_dcy_param ! routine used here and called by warm-layer parameterization (sbcblk_skin_coare*) 31 38 32 39 !!---------------------------------------------------------------------- 33 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 34 !! $Id$ 41 !! $Id$ 35 42 !! Software governed by the CeCILL license (see ./LICENSE) 36 43 !!---------------------------------------------------------------------- 37 44 CONTAINS 38 45 39 40 41 42 43 44 & rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )45 46 47 48 46 INTEGER FUNCTION sbc_dcy_alloc() 47 !!---------------------------------------------------------------------- 48 !! *** FUNCTION sbc_dcy_alloc *** 49 !!---------------------------------------------------------------------- 50 ALLOCATE( raa (jpi,jpj) , rbb (jpi,jpj) , rcc (jpi,jpj) , rab (jpi,jpj) , & 51 & rtmd(jpi,jpj) , rdawn_dcy(jpi,jpj) , rdusk_dcy(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 52 ! 53 CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 54 IF( sbc_dcy_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 55 END FUNCTION sbc_dcy_alloc 49 56 50 57 … … 60 67 !! 61 68 !! reference : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007 62 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 69 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 63 70 !! Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 64 71 !!---------------------------------------------------------------------- 65 72 LOGICAL , OPTIONAL , INTENT(in) :: l_mask ! use the routine for night mask computation 66 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqsrin ! input daily QSR flux 73 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqsrin ! input daily QSR flux 67 74 REAL(wp), DIMENSION(jpi,jpj) :: zqsrout ! output QSR flux with diurnal cycle 68 75 !! 69 76 INTEGER :: ji, jj ! dummy loop indices 70 77 INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 71 REAL(wp) :: ztwopi, zinvtwopi, zconvrad72 78 REAL(wp) :: zlo, zup, zlousd, zupusd 73 REAL(wp) :: zdsws, zdecrad, ztx, zsin, zcos 74 REAL(wp) :: ztmp, ztmp1, ztmp2, ztest 79 REAL(wp) :: ztmp, ztmp1, ztmp2 75 80 REAL(wp) :: ztmpm, ztmpm1, ztmpm2 76 !---------------------------statement functions------------------------77 REAL(wp) :: fintegral, pt1, pt2, paaa, pbbb, pccc ! dummy statement function arguments78 fintegral( pt1, pt2, paaa, pbbb, pccc ) = &79 & paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2) &80 & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)81 81 !!--------------------------------------------------------------------- 82 82 ! 83 83 ! Initialization 84 84 ! -------------- 85 ztwopi = 2._wp * rpi86 zinvtwopi = 1._wp / ztwopi87 zconvrad = ztwopi / 360._wp88 89 85 ! When are we during the day (from 0 to 1) 90 86 zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday 91 87 zup = zlo + ( REAL(nn_fsbc, wp) * rdt ) / rday 92 ! 93 IF( nday_qsr == -1 ) THEN ! first time step only 88 ! 89 IF( nday_qsr == -1 ) THEN ! first time step only 94 90 IF(lwp) THEN 95 91 WRITE(numout,*) … … 98 94 WRITE(numout,*) 99 95 ENDIF 96 ENDIF 97 98 ! Setting parameters for each new day: 99 CALL sbc_dcy_param() 100 101 !CALL iom_put( "rdusk_dcy", rdusk_dcy(:,:)*tmask(:,:,1) ) !LB 102 !CALL iom_put( "rdawn_dcy", rdawn_dcy(:,:)*tmask(:,:,1) ) !LB 103 !CALL iom_put( "rscal_dcy", rscal(:,:)*tmask(:,:,1) ) !LB 104 105 106 ! 3. update qsr with the diurnal cycle 107 ! ------------------------------------ 108 109 imask_night(:,:) = 0 110 DO jj = 1, jpj 111 DO ji = 1, jpi 112 ztmpm = 0._wp 113 IF( ABS(rab(ji,jj)) < 1. ) THEN ! day duration is less than 24h 114 ! 115 IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN ! day time in one part 116 zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 117 zlousd = MIN(zlousd, zup) 118 zupusd = MIN(zup, rdusk_dcy(ji,jj)) 119 zupusd = MAX(zupusd, zlo) 120 ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 121 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 122 ztmpm = zupusd - zlousd 123 IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 124 ! 125 ELSE ! day time in two parts 126 zlousd = MIN(zlo, rdusk_dcy(ji,jj)) 127 zupusd = MIN(zup, rdusk_dcy(ji,jj)) 128 ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 129 ztmpm1=zupusd-zlousd 130 zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 131 zupusd = MAX(zup, rdawn_dcy(ji,jj)) 132 ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 133 ztmpm2 =zupusd-zlousd 134 ztmp = ztmp1 + ztmp2 135 ztmpm = ztmpm1 + ztmpm2 136 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 137 IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1 138 ENDIF 139 ELSE ! 24h light or 24h night 140 ! 141 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 142 ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 143 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 144 imask_night(ji,jj) = 0 145 ! 146 ELSE ! No day 147 zqsrout(ji,jj) = 0.0_wp 148 imask_night(ji,jj) = 1 149 ENDIF 150 ENDIF 151 END DO 152 END DO 153 ! 154 IF( PRESENT(l_mask) .AND. l_mask ) THEN 155 zqsrout(:,:) = float(imask_night(:,:)) 156 ENDIF 157 ! 158 END FUNCTION sbc_dcy 159 160 161 SUBROUTINE sbc_dcy_param( ) 162 !! 163 INTEGER :: ji, jj ! dummy loop indices 164 !INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 165 REAL(wp) :: zdsws, zdecrad, ztx, zsin, zcos 166 REAL(wp) :: ztmp, ztest 167 !---------------------------statement functions------------------------ 168 ! 169 IF( nday_qsr == -1 ) THEN ! first time step only 100 170 ! allocate sbcdcy arrays 101 171 IF( sbc_dcy_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' ) 102 172 ! Compute rcc needed to compute the time integral of the diurnal cycle 103 rcc(:,:) = zconvrad * glamt(:,:) - rpi173 rcc(:,:) = rad * glamt(:,:) - rpi 104 174 ! time of midday 105 175 rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp … … 107 177 ENDIF 108 178 109 ! If this is a new day, we have to update the dawn, dusk and scaling function 179 ! If this is a new day, we have to update the dawn, dusk and scaling function 110 180 !---------------------- 111 112 ! 2.1 dawn and dusk 113 114 ! nday is the number of days since the beginning of the current month 115 IF( nday_qsr /= nday ) THEN 181 182 ! 2.1 dawn and dusk 183 184 ! nday is the number of days since the beginning of the current month 185 IF( nday_qsr /= nday ) THEN 116 186 ! save the day of the year and the daily mean of qsr 117 nday_qsr = nday 118 ! number of days since the previous winter solstice (supposed to be always 21 December) 187 nday_qsr = nday 188 ! number of days since the previous winter solstice (supposed to be always 21 December) 119 189 zdsws = REAL(11 + nday_year, wp) 120 190 ! declination of the earths orbit 121 zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )191 zdecrad = (-23.5_wp * rad) * COS( zdsws * 2._wp*rpi / REAL(nyear_len(1),wp) ) 122 192 ! Compute A and B needed to compute the time integral of the diurnal cycle 123 193 … … 125 195 DO jj = 1, jpj 126 196 DO ji = 1, jpi 127 ztmp = zconvrad * gphit(ji,jj)197 ztmp = rad * gphit(ji,jj) 128 198 raa(ji,jj) = SIN( ztmp ) * zsin 129 199 rbb(ji,jj) = COS( ztmp ) * zcos 130 END DO 131 END DO 200 END DO 201 END DO 132 202 ! Compute the time of dawn and dusk 133 203 134 ! rab to test if the day time is equal to 0, less than 24h of full day 204 ! rab to test if the day time is equal to 0, less than 24h of full day 135 205 rab(:,:) = -raa(:,:) / rbb(:,:) 136 206 DO jj = 1, jpj 137 207 DO ji = 1, jpi 138 IF 139 ! When is it night?140 ztx = zinvtwopi* (ACOS(rab(ji,jj)) - rcc(ji,jj))141 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )142 ! is it dawn or dusk?143 IF 144 rdawn (ji,jj) = ztx145 rdusk (ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )208 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 209 ! When is it night? 210 ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 211 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 212 ! is it dawn or dusk? 213 IF( ztest > 0._wp ) THEN 214 rdawn_dcy(ji,jj) = ztx 215 rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) 146 216 ELSE 147 rdusk (ji,jj) = ztx148 rdawn (ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )217 rdusk_dcy(ji,jj) = ztx 218 rdawn_dcy(ji,jj) = rtmd(ji,jj) - ( rdusk_dcy(ji,jj) - rtmd(ji,jj) ) 149 219 ENDIF 150 220 ELSE 151 rdawn (ji,jj) = rtmd(ji,jj) + 0.5_wp152 rdusk (ji,jj) = rdawn(ji,jj)153 ENDIF 154 END DO155 END DO 156 rdawn (:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )157 rdusk (:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )221 rdawn_dcy(ji,jj) = rtmd(ji,jj) + 0.5_wp 222 rdusk_dcy(ji,jj) = rdawn_dcy(ji,jj) 223 ENDIF 224 END DO 225 END DO 226 rdawn_dcy(:,:) = MOD( (rdawn_dcy(:,:) + 1._wp), 1._wp ) 227 rdusk_dcy(:,:) = MOD( (rdusk_dcy(:,:) + 1._wp), 1._wp ) 158 228 ! 2.2 Compute the scaling function: 159 229 ! S* = the inverse of the time integral of the diurnal cycle from dawn to dusk … … 162 232 DO jj = 1, jpj 163 233 DO ji = 1, jpi 164 IF 234 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 165 235 rscal(ji,jj) = 0.0_wp 166 IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part167 IF( (rdusk (ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN168 rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))169 rscal(ji,jj) = 1._wp / rscal(ji,jj)236 IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN ! day time in one part 237 IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 238 rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 239 rscal(ji,jj) = 1._wp / rscal(ji,jj) 170 240 ENDIF 171 241 ELSE ! day time in two parts 172 IF( (rdusk (ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN173 rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) &174 & + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))175 rscal(ji,jj) = 1. / rscal(ji,jj)242 IF( (rdusk_dcy(ji,jj) + (1._wp - rdawn_dcy(ji,jj)) ) .ge. 0.001_wp ) THEN 243 rscal(ji,jj) = fintegral(0._wp, rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) & 244 & + fintegral(rdawn_dcy(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 245 rscal(ji,jj) = 1. / rscal(ji,jj) 176 246 ENDIF 177 247 ENDIF 178 248 ELSE 179 IF 180 rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 249 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 250 rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 181 251 rscal(ji,jj) = 1._wp / rscal(ji,jj) 182 252 ELSE ! No day … … 184 254 ENDIF 185 255 ENDIF 186 END DO 187 END DO 256 END DO 257 END DO 188 258 ! 189 259 ztmp = rday / ( rdt * REAL(nn_fsbc, wp) ) 190 260 rscal(:,:) = rscal(:,:) * ztmp 191 261 ! 192 ENDIF 193 ! 3. update qsr with the diurnal cycle 194 ! ------------------------------------ 195 196 imask_night(:,:) = 0 197 DO jj = 1, jpj 198 DO ji = 1, jpi 199 ztmpm = 0._wp 200 IF( ABS(rab(ji,jj)) < 1. ) THEN ! day duration is less than 24h 201 ! 202 IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part 203 zlousd = MAX(zlo, rdawn(ji,jj)) 204 zlousd = MIN(zlousd, zup) 205 zupusd = MIN(zup, rdusk(ji,jj)) 206 zupusd = MAX(zupusd, zlo) 207 ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 208 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 209 ztmpm = zupusd - zlousd 210 IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 211 ! 212 ELSE ! day time in two parts 213 zlousd = MIN(zlo, rdusk(ji,jj)) 214 zupusd = MIN(zup, rdusk(ji,jj)) 215 ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 216 ztmpm1=zupusd-zlousd 217 zlousd = MAX(zlo, rdawn(ji,jj)) 218 zupusd = MAX(zup, rdawn(ji,jj)) 219 ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 220 ztmpm2 =zupusd-zlousd 221 ztmp = ztmp1 + ztmp2 222 ztmpm = ztmpm1 + ztmpm2 223 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 224 IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 225 ENDIF 226 ELSE ! 24h light or 24h night 227 ! 228 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 229 ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 230 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 231 imask_night(ji,jj) = 0 232 ! 233 ELSE ! No day 234 zqsrout(ji,jj) = 0.0_wp 235 imask_night(ji,jj) = 1 236 ENDIF 237 ENDIF 238 END DO 239 END DO 240 ! 241 IF( PRESENT(l_mask) .AND. l_mask ) THEN 242 zqsrout(:,:) = float(imask_night(:,:)) 243 ENDIF 244 ! 245 END FUNCTION sbc_dcy 262 ENDIF !IF( nday_qsr /= nday ) 263 ! 264 END SUBROUTINE sbc_dcy_param 265 266 267 FUNCTION fintegral( pt1, pt2, paaa, pbbb, pccc ) 268 REAL(wp), INTENT(in) :: pt1, pt2, paaa, pbbb, pccc 269 REAL(wp) :: fintegral 270 fintegral = paaa * pt2 + 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt2) & 271 & - paaa * pt1 - 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt1) 272 END FUNCTION fintegral 246 273 247 274 !!====================================================================== -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcmod.F90
r11874 r12015 128 128 ENDIF 129 129 #else 130 IF( lk_si3 ) nn_ice = 2130 !IF( lk_si3 ) nn_ice = 2 131 131 IF( lk_cice ) nn_ice = 3 132 132 #endif … … 261 261 262 262 ! ! Choice of the Surface Boudary Condition (set nsbc) 263 nday_qsr = -1 ! allow initialization at the 1st call !LB: now warm-layer of COARE* calls "sbc_dcy_param" of sbcdcy.F90! 263 264 IF( ln_dm2dc ) THEN !* daily mean to diurnal cycle 264 nday_qsr = -1 ! allow initialization at the 1st call265 !LB:nday_qsr = -1 ! allow initialization at the 1st call 265 266 IF( .NOT.( ln_flx .OR. ln_blk .OR. ln_abl ) .AND. nn_components /= jp_iam_opa ) & 266 267 & CALL ctl_stop( 'qsr diurnal cycle from daily values requires flux, bulk or abl formulation' ) … … 370 371 CALL iom_set_rstw_var_active('sfx_b') 371 372 ENDIF 372 ! 373 373 374 END SUBROUTINE sbc_init 374 375 … … 409 410 emp_b (:,:) = emp (:,:) 410 411 sfx_b (:,:) = sfx (:,:) 411 IF 412 IF( ln_rnf ) THEN 412 413 rnf_b (:,: ) = rnf (:,: ) 413 414 rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) … … 451 452 IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! forced-coupled mixed formulation after forcing 452 453 ! 453 IF 454 IF( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( ) ! Wind stress provided by waves 454 455 ! 455 456 ! !== Misc. Options ==! … … 485 486 !!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 486 487 !!$ CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) 487 IF 488 IF( ll_wd ) THEN ! If near WAD point limit the flux for now 488 489 zthscl = atanh(rn_wd_sbcfra) ! taper frac default is .999 489 490 zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1 ! do this calc of water
Note: See TracChangeset
for help on using the changeset viewer.
