pbl_surface Subroutine

public subroutine pbl_surface(dtime, date0, itap, jour, debut, lafin, rlon, rlat, rugoro, rmu0, lwdown_m, pphi, cldt, rain_f, snow_f, bs_f, solsw_m, solswfdiff_m, sollw_m, gustiness, t, q, qbs, u, v, wake_dlt, wake_dlq, wake_cstar, wake_s, pplay, paprs, pctsrf, ts, SFRWL, alb_dir, alb_dif, ustar, u10m, v10m, wstar, cdragh, cdragm, zu1, zv1, beta, alb_dir_m, alb_dif_m, zxsens, zxevap, zxsnowerosion, icesub_lic, alb3_lic, runoff, snowhgt, qsnow, to_ice, sissnow, zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, d_t, d_q, d_qbs, d_u, d_v, d_t_diss, d_t_w, d_q_w, d_t_x, d_q_x, zxsens_x, zxfluxlat_x, zxsens_w, zxfluxlat_w, delta_tsurf, wake_dens, cdragh_x, cdragh_w, cdragm_x, cdragm_w, kh, kh_x, kh_w, zcoefh, zcoefm, slab_wfbils, qsol, zq2m, s_pblh, s_plcl, s_pblh_x, s_plcl_x, s_pblh_w, s_plcl_w, s_capCL, s_oliqCL, s_cteiCL, s_pblT, s_therm, s_trmb1, s_trmb2, s_trmb3, zustar, zu10m, zv10m, fder_print, zxqsurf, delta_qsurf, rh2m, zxfluxu, zxfluxv, z0m, z0h, agesno, sollw, solsw, d_ts, evap, fluxlat, t2m, wfbils, wfevap, flux_t, flux_u, flux_v, dflux_t, dflux_q, zxsnow, zxfluxt, zxfluxq, zxfluxqbs, q2m, flux_q, flux_qbs, tke_x, eps_x, wake_dltke, treedrg, hice, tice, bilg_cumul, fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, dtice_melt, dtice_snow2sic, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, cdragm_tersrf, cdragh_tersrf, swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf)

! nrlmd+jyg le 02/05/2011 et le 20/02/2012 t_x, q_x, t_w, q_w, & ! jyg ! nrlmd+jyg le 02/05/2011 et le 20/02/2012 d_wake_dlt,d_wake_dlq, & ! ! nrlmd le 13/06/2011 ! ! ! jyg le 08/02/2012 ! zxfluxt, zxfluxq, q2m, flux_q, tke, & jyg ! nrlmd+jyg le 02/05/2011 et le 20/02/2012 tke_x, tke_w & !

Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 Objet: interface de "couche limite" (diffusion verticale)

AA REM: AA----- AA Tout ce qui a trait au traceurs est dans phytrac maintenant AA pour l'instant le calcul de la couche limite pour les traceurs AA se fait avec cltrac et ne tient pas compte de la differentiation AA des sous-fraction de sol. AA REM bis : AA---------- AA Pour pouvoir extraire les coefficient d'echanges et le vent AA dans la premiere couche, 3 champs supplementaires ont ete crees AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir AA si les informations des subsurfaces doivent etre prises en compte AA il faudra sortir ces memes champs en leur ajoutant une dimension, AA c'est a dire nbsrf (nbre de subsurface).

Arguments:

dtime----input-R- interval du temps (secondes) itap-----input-I- numero du pas de temps date0----input-R- jour initial t--------input-R- temperature (K) q--------input-R- vapeur d'eau (kg/kg) u--------input-R- vitesse u v--------input-R- vitesse v wake_dlt-input-R- temperatre difference between (w) and (x) (K) wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg) wake_cstar-input-R- wake gust front speed (m/s) wake_s---input-R- wake fractionnal area ts-------input-R- temperature du sol (en Kelvin) paprs----input-R- pression a intercouche (Pa) pplay----input-R- pression au milieu de couche (Pa) rlat-----input-R- latitude en degree z0m, z0h ----input-R- longeur de rugosite (en m) Martin cldt-----input-R- total cloud fraction Martin GG pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol) GG

d_t------output-R- le changement pour "t" d_q------output-R- le changement pour "q" d_u------output-R- le changement pour "u" d_v------output-R- le changement pour "v" d_ts-----output-R- le changement pour "ts" flux_t---output-R- flux de chaleur sensible (CpT) J/m2/s (W/m2) (orientation positive vers le bas) tke_x---input/output-R- tke in the (x) region (kg/m2/s) wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m2/s) flux_q---output-R- flux de vapeur d'eau (kg/m2/s) flux_u---output-R- tension du vent X: (kg m/s)/(m2 s) ou Pascal flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal dflux_t--output-R- derive du flux sensible dflux_q--output-R- derive du flux latent zu1------output-R- le vent dans la premiere couche zv1------output-R- le vent dans la premiere couche trmb1----output-R- deep_cape trmb2----output-R- inhibition trmb3----output-R- Point Omega cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL plcl-----output-R- Niveau de condensation pblh-----output-R- HCL pblT-----output-R- T au nveau HCL treedrg--output-R- tree drag (m)
qsurf_tersrf--output-R- surface specific humidity of continental sub-surfaces cdragm_tersrf--output-R- momentum drag coefficient of continental sub-surfaces cdragh_tersrf--output-R- heat drag coefficient of continental sub-surfaces tsurf_new_tersrf--output-R- surface temperature of continental sub-surfaces swnet_tersrf--output-R- net shortwave radiation of continental sub-surfaces lwnet_tersrf--output-R- net longwave radiation of continental sub-surfaces fluxsens_tersrf--output-R- sensible heat flux of continental sub-surfaces fluxlat_tersrf--output-R- latent heat flux of continental sub-surfaces

jyg iflag_split = mod(iflag_pbl_split,2) jyg iflag_split = mod(iflag_pbl_split,10)

1) Initialisation and validation tests Only done first time entering this subroutine

Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket instead of ORCHIDEE)

2) Initialization to zero

2a) Initialization of all argument variables with INTENT(OUT)

tke(:,:,is_ave)=0.

jyg ! jyg le 23/02/2013 ! ! ! jyg le 10/02/2012

cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
zv1 = 0.0     ; yqsurf = 0.0

d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0

! nrlmd+jyg le 02/05/2011 et le 20/02/2012 d_t_w=0. ; d_q_w=0.
d_t_x=0. ; d_q_x=0. d_wake_dlt=0. ; d_wake_dlq=0. ! ! nrlmd le 13/06/2011 !

jyg jyg !

!

3) - Calculate pressure thickness of each layer - Calculate the wind at first layer - Mean calculations of albedo - Calculate net radiance at sub-surface

Test for rugos........ from physiq.. A la fin plutot???

al1

4) Loop over different surfaces

Only points containing a fraction of the sub surface will be treated.

<<<<<<<<<<<<< (nsrf == is_oce) (iflag_split_ref == 3)

! jyg le 19/08/2012

!

5) Compress variables

jyg ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 ) jyg ! nrlmd le 13/06/2011 ! ! jyg le 07/02/2012 et le 10/04/2013 ytke(j,k) = tke(i,k,nsrf) jyg GG

! nrlmd le 02/05/2011 !

! nrlmd le 02/05/2011 ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)wake_dltke(i,k,nsrf) ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))wake_dltke(i,k,nsrf) ywake_dltke(j,k) = wake_dltke(i,k,nsrf) ytke(j,k) = tke(i,k,nsrf)

jyg ! ! jyg le 07/02/2012 ! ! nrlmd le 13/06/2011 ! !

6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.

! jyg le 07/02/2012 ! ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012

! AM heterogeneous continental subsurfaces ! !

6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.

! jyg le 07/02/2012 ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 !jyg le 10/04/2013 En attendant de traiter le transport des traceurs dans les poches froides, formule arbitraire pour ycoefh et ycoefm

8) "La descente" - "The downhill"

climb_hq_down and climb_wind_down calculate the coefficients Ccoef_X et Dcoef_X for X=[H, Q, U, V]. Only the coefficients at surface for H and Q are returned.

! jyg le 07/02/2012 ! ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 ! jyg le 09/05/2011 ! ! nrlmd le 02/05/2011 ! ! ! nrlmd le 02/05/2011 ! ! ! !

! jyg le 07/02/2012 ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 ! jyg le 09/05/2011 ! ! nrlmd le 02/05/2011 !

! nrlmd le 02/05/2011 ! !
!

9) Small calculations

! nrlmd le 02/05/2011 -----------------------On raccorde les 2 colonnes dans la couche 1

! jyg le 07/02/2012 ! jyg le 01/02/2017

Cdragq computed from cdrag The difference comes only from a factor (f_z0qh_oce) on z0, so that it can be computed inside wx_pbl0_merge More complicated appraches may require the propagation through pbl_surface of an independant cdragq variable.

   ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
        log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
   ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
        log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))

 Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag

Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', & z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon) !

Calulate t2m and q2m for the case of calculation at land grid points t2m and q2m are needed as input to ORCHIDEE

10) Switch according to current surface It is necessary to start with the continental surfaces because the ocean needs their run-off.

jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&

      alb3_lic(:)=0.

jyg jyg ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),& jyg yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&

11) - Calcul the increment of surface temperature

12) "La remontee" - "The uphill"

The fluxes (y_flux_X) and tendancy (y_d_X) are calculated for X=H, Q, U and V, for all vertical levels.

! ! jyg le 10/04/2013 et EV 10/2020

! Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg) IF (iflag_split .eq.0) THEN ENDIF ! (iflag_split .eq.0)

!! HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &

jyg jyg!! A reprendre apres reflexion =============================================== jyg!! jyg!! DO j=1,knon jyg!!!!! nrlmd le 13/06/2011 jyg!! jyg!!!----Diffusion dans le sol dans le cas continental seulement jyg!! IF (nsrf.eq.is_ter) THEN jyg!!!----Calcul du coefficient delta_coeff jyg!! tau_eq(j)=(ywake_s(j)/2.)(1./max(wake_cstar(j),0.01))sqrt(0.4/(3.14max(wake_dens(j),8e-12))) jyg!! jyg!!! delta_coef(j)=dtime/(inertiasqrt(tau_eq(j))) jyg!! delta_coef(j)=facteursqrt(tau_eq(j))/inertia jyg!!! delta_coef(j)=0. jyg!! ELSE jyg!! delta_coef(j)=0. jyg!! ENDIF jyg!! jyg!!!----Calcul de delta_tsurf jyg!! y_delta_tsurf(j)=delta_coef(j)y_delta_flux_t1(j) jyg!! jyg!!!----Si il n'y a pas des poches... jyg!! IF (wake_cstar(j).le.0.01) THEN jyg!! y_delta_tsurf(j)=0. jyg!! y_delta_flux_t1(j)=0. jyg!! ENDIF jyg!! jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle) jyg!!!!!!! jyg le 23/02/2012 jyg!!!!!!! jyg!!!! ybeta(j)=y_flux_q1(j) / & jyg!!!! & (Kech_h(j)(yq(j,1)-yqsatsurf(j))) jyg!!!!!! ybeta(j)=-1.yevap(j) / & jyg!!!!!! & (ywake_s(j)Kech_h_w(j)(yq_w(j,1)-yqsatsurf_w(j))+(1.-ywake_s(j))Kech_h_x(j)(yq_x(j,1)-yqsatsurf_x(j))) jyg!!!!!!! fin jyg jyg!!!!! jyg!! jyg!! ENDDO jyg!! jyg!!!!! fin nrlmd le 13/06/2011 jyg!! ! jyg le 07/02/2012 !

! jyg le 07/02/2012 ! ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 ! jyg le 07/02/2012 ! ! nrlmd le 02/05/2011 !

! nrlmd le 02/05/2011 ! ! !

! jyg le 07/02/2012 ! ! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012 ! jyg le 07/02/2012 ! ! nrlmd le 02/05/2011 !

! nrlmd le 02/05/2011 ! ! !

    DO j = 1, knon
      y_dflux_t(j) = y_dflux_t(j) * ypct(j)
      y_dflux_q(j) = y_dflux_q(j) * ypct(j)
    ENDDO

13) Transform variables for output format : - Decompress - Multiply with pourcentage of current surface - Cumulate in global variable

! jyg le 07/02/2012 ! ! jyg le 07/02/2012 ! ! nrlmd le 02/05/2011 ! ! nrlmd le 13/06/2011 jyg20170131 delta_tsurf(i,nsrf)=y_delta_tsurf(j)ypct(j) jyg20210118 delta_tsurf(i,nsrf)=y_delta_tsurf(j) ! !
! ! nrlmd le 02/05/2011 jyg le 20/02/2011 tke_x(:,:,nsrf)=0. tke_w(:,:,nsrf)=0. jyg le 20/02/2011 DO k = 1, klev+1 DO j = 1, knon i = ni(j) wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k) tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)wake_dltke(i,k,nsrf) ENDDO ENDDO jyg le 20/02/2011 DO k = 1, klev+1 DO j = 1, knon i = ni(j) tke(i,k,nsrf)=(1.-ywake_s(j))tke_x(i,k,nsrf)+ywake_s(j)tke_w(i,k,nsrf) ENDDO ENDDO ! tke(i,k,nsrf) = ytke(j,k) tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)ypct(j) jyg tke(i,k,nsrf) = ytke_x(j,k) + ywake_s(j)wake_dltke(i,k,nsrf) tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j) jyg ! ftsoil(:,:,nsrf) = 0. jyg ! jyg le 07/02/2012 ! ! nrlmd+jyg le 02/05/2011 et le 20/02/2012 d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k) d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k) ! !

14) Calculate the temperature and relative humidity at 2m and the wind at 10m Call HBTM

!

! jyg le 07/02/2012 tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j) !
! ! jyg le 07/02/2012 ! ! ! jyg le 07/02/2012 ! !

! jyg le 07/02/2012 !
! ! jyg le 07/02/2012 !
!

! jyg le 07/02/2012 ! ! !

15) End of loop over different surfaces

16) Calculate the mean value over all sub-surfaces for some variables

! jyg le 07/02/2012 ! ! nrlmd & jyg les 02/05/2011, 05/02/2012

! !

!

! jyg le 07/02/2012 ! ! nrlmd le 02/05/2011 !

al1 ! jyg le 07/02/2012 ! nrlmd le 02/05/2011 ! ! jyg le 08/02/2012 Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ; pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ; pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ; pour les autres variables, on sort les valeurs de la region (x). ! !

$! $! If a sub-surface does not exsist for a grid point, the mean value for all $! sub-surfaces is distributed. $! $ DO nsrf = 1, nbsrf $ DO i = 1, klon $ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN $ ts(i,nsrf) = zxtsol(i) $ t2m(i,nsrf) = zt2m(i) $ q2m(i,nsrf) = zq2m(i) $ u10m(i,nsrf) = zu10m(i) $ v10m(i,nsrf) = zv10m(i) $ $! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour $ pblh(i,nsrf) = s_pblh(i) $ plcl(i,nsrf) = s_plcl(i) $ capCL(i,nsrf) = s_capCL(i) $ oliqCL(i,nsrf) = s_oliqCL(i) $ cteiCL(i,nsrf) = s_cteiCL(i) $ pblT(i,nsrf) = s_pblT(i) $ therm(i,nsrf) = s_therm(i) $ trmb1(i,nsrf) = s_trmb1(i) $ trmb2(i,nsrf) = s_trmb2(i) $ trmb3(i,nsrf) = s_trmb3(i) $ ENDIF $ ENDDO $ ENDDO

Arguments