! ! $Header$ ! MODULE pbl_surface_mod ! ! Planetary Boundary Layer and Surface module ! ! This module manage the calculation of turbulent diffusion in the boundary layer ! and all interactions towards the differents sub-surfaces. ! ! USE dimphy USE mod_phys_lmdz_para, ONLY : mpi_size USE ioipsl USE surface_data, ONLY : ocean, ok_veget, debug_surf, newwind USE surf_land_mod, ONLY : surf_land USE surf_landice_mod, ONLY : surf_landice USE surf_ocean_mod, ONLY : surf_ocean USE surf_seaice_mod, ONLY : surf_seaice USE cpl_mod, ONLY : gath2cpl USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up USE coef_diff_turb_mod, ONLY : coef_diff_turb !jg+ temporary USE mod_clvent, ONLY : clvent, save_flux !jg- IMPLICIT NONE ! Declaration of variables saved in restart file REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol !$OMP THREADPRIVATE(qsol) REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder !$OMP THREADPRIVATE(fder) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow !$OMP THREADPRIVATE(snow) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf !$OMP THREADPRIVATE(qsurf) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap !$OMP THREADPRIVATE(evap) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos !$OMP THREADPRIVATE(rugos) REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno !$OMP THREADPRIVATE(agesno) REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil !$OMP THREADPRIVATE(ftsoil) CONTAINS ! !**************************************************************************************** ! SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,& evap_rst, rugos_rst, agesno_rst, ftsoil_rst) ! This routine should be called after the restart file has been read. ! This routine initialize the restart variables and does some validation tests ! for the index of the different surfaces and tests the choice of type of ocean. INCLUDE "indicesol.inc" INCLUDE "dimsoil.h" INCLUDE "iniprint.h" ! Input variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst REAL, DIMENSION(klon), INTENT(IN) :: fder_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst ! Local variables !**************************************************************************************** INTEGER :: ierr CHARACTER(len=80) :: abort_message CHARACTER(len = 20) :: modname = 'pbl_surface_init' !**************************************************************************************** ! Allocate and initialize module variables with fields read from restart file. ! !**************************************************************************************** ALLOCATE(qsol(klon), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(fder(klon), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(snow(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(qsurf(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(evap(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(rugos(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(agesno(klon,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr) IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1) qsol(:) = qsol_rst(:) fder(:) = fder_rst(:) snow(:,:) = snow_rst(:,:) qsurf(:,:) = qsurf_rst(:,:) evap(:,:) = evap_rst(:,:) rugos(:,:) = rugos_rst(:,:) agesno(:,:) = agesno_rst(:,:) ftsoil(:,:,:) = ftsoil_rst(:,:,:) !**************************************************************************************** ! Test for sub-surface indices ! !**************************************************************************************** IF (is_ter /= 1) THEN WRITE(lunout,*)" *** Warning ***" WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter WRITE(lunout,*)"or on doit commencer par les surfaces continentales" abort_message="voir ci-dessus" CALL abort_gcm(modname,abort_message,1) ENDIF IF ( is_oce > is_sic ) THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)' Pour des raisons de sequencement dans le code' WRITE(lunout,*)' l''ocean doit etre traite avant la banquise' WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic abort_message='voir ci-dessus' CALL abort_gcm(modname,abort_message,1) ENDIF IF ( is_lic > is_sic ) THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)' Pour des raisons de sequencement dans le code' WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer' WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic abort_message='voir ci-dessus' CALL abort_gcm(modname,abort_message,1) ENDIF !**************************************************************************************** ! Validation of ocean mode ! !**************************************************************************************** IF (ocean /= 'slab ' .AND. ocean /= 'force ' .AND. ocean /= 'couple') THEN WRITE(lunout,*)' *** Warning ***' WRITE(lunout,*)'Option couplage pour l''ocean = ', ocean abort_message='option pour l''ocean non valable' CALL abort_gcm(modname,abort_message,1) ENDIF !**************************************************************************************** ! Test of coherence between variable ok_veget and cpp key CPP_VEGET ! !**************************************************************************************** IF (ok_veget) THEN #ifndef CPP_VEGET abort_message='Pb de coherence: ok_veget = .true. mais CPP_VEGET = .false.' CALL abort_gcm(modname,abort_message,1) #endif ENDIF END SUBROUTINE pbl_surface_init ! !**************************************************************************************** ! SUBROUTINE pbl_surface( & dtime, date0, itap, jour, & debut, lafin, & rlon, rlat, rugoro, rmu0, & rain_f, snow_f, solsw_m, sollw_m, & t, q, u, v, & pplay, paprs, pctsrf, & ts, albe, alblw, u10m, v10m, & sollwdown, cdragh, cdragm, zu1, zv1, & albsol, albsollw, zxsens, zxevap, & zxtsol, zxfluxlat, zt2m, qsat2m, & d_t, d_q, d_u, d_v, & zcoefh, pctsrf_new, & qsol_d, zq2m, s_pblh, s_plcl, & s_capCL, s_oliqCL, s_cteiCL, s_pblT, & s_therm, s_trmb1, s_trmb2, s_trmb3, & zxrugs, zu10m, zv10m, fder_print, & zxqsurf, rh2m, zxfluxu, zxfluxv, & rugos_d, agesno_d, sollw, solsw, & d_ts, evap_d, fluxlat, t2m, & wfbils, wfbilo, flux_t, flux_u, flux_v,& dflux_t, dflux_q, zxsnow, & zxfluxt, zxfluxq, q2m, flux_q ) !**************************************************************************************** ! 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 ! 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 ! rugos----input-R- longeur de rugosite (en m) ! ! 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/m**2/s (W/m**2) ! (orientation positive vers le bas) ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s) ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 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 ! INCLUDE "indicesol.inc" INCLUDE "dimsoil.h" INCLUDE "YOMCST.inc" INCLUDE "iniprint.h" INCLUDE "FCTTRE.inc" INCLUDE "clesphys.inc" INCLUDE "compbl.h" INCLUDE "dimensions90.h" INCLUDE "YOETHF.inc" INCLUDE "temps.inc" INCLUDE "control.inc" ! Input variables !**************************************************************************************** REAL, INTENT(IN) :: dtime REAL, INTENT(IN) :: date0 INTEGER, INTENT(IN) :: itap INTEGER, INTENT(IN) :: jour ! jour de l'annee en cours LOGICAL, INTENT(IN) :: debut, lafin REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat REAL, DIMENSION(klon), INTENT(IN) :: rugoro REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosinus de l'angle solaire zenithal REAL, DIMENSION(klon), INTENT(IN) :: rain_f, snow_f REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! mean value REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! mean value REAL, DIMENSION(klon,klev), INTENT(IN) :: t, q REAL, DIMENSION(klon,klev), INTENT(IN) :: u, v REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! Input/Output variables !**************************************************************************************** REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: albe REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alblw REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m, v10m ! Output variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: sollwdown REAL, DIMENSION(klon), INTENT(OUT) :: cdragh, cdragm REAL, DIMENSION(klon), INTENT(OUT) :: zu1 REAL, DIMENSION(klon), INTENT(OUT) :: zv1 REAL, DIMENSION(klon), INTENT(OUT) :: albsol REAL, DIMENSION(klon), INTENT(OUT) :: albsollw REAL, DIMENSION(klon), INTENT(OUT) :: zxsens, zxevap REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat REAL, DIMENSION(klon), INTENT(OUT) :: zt2m REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t, d_q REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u, d_v REAL, DIMENSION(klon, klev), INTENT(OUT) :: zcoefh REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: pctsrf_new ! Output only for diagnostics REAL, DIMENSION(klon), INTENT(OUT) :: qsol_d REAL, DIMENSION(klon), INTENT(OUT) :: zq2m REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl REAL, DIMENSION(klon), INTENT(OUT) :: s_capCL REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqCL REAL, DIMENSION(klon), INTENT(OUT) :: s_cteiCL REAL, DIMENSION(klon), INTENT(OUT) :: s_pblT REAL, DIMENSION(klon), INTENT(OUT) :: s_therm REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 REAL, DIMENSION(klon), INTENT(OUT) :: zxrugs REAL, DIMENSION(klon), INTENT(OUT) :: zu10m REAL, DIMENSION(klon), INTENT(OUT) :: zv10m REAL, DIMENSION(klon), INTENT(OUT) :: fder_print REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf REAL, DIMENSION(klon), INTENT(OUT) :: rh2m REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu, zxfluxv REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_d REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_d REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw, solsw REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_d REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils, wfbilo REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u, flux_v ! Output not needed REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t, dflux_q REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt, zxfluxq REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! Local variables with attribute SAVE !**************************************************************************************** INTEGER, SAVE :: nhoridbg, nidbg !$OMP THREADPRIVATE(nhoridbg, nidbg) LOGICAL, SAVE :: debugindex=.FALSE. !$OMP THREADPRIVATE(debugindex) LOGICAL, SAVE :: first_call=.TRUE. !$OMP THREADPRIVATE(first_call) CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf !$OMP THREADPRIVATE(cl_surf) ! Other local variables !**************************************************************************************** INTEGER :: i, k, nsrf INTEGER :: knon, j INTEGER :: idayref INTEGER , DIMENSION(klon) :: ni REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola REAL :: amn, amx REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere REAL, DIMENSION(klon) :: yts, yrugos, ypct, yz0_new REAL, DIMENSION(klon) :: yalb REAL, DIMENSION(klon) :: yalblw REAL, DIMENSION(klon) :: yu1, yv1 REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol REAL, DIMENSION(klon) :: yrain_f, ysnow_f REAL, DIMENSION(klon) :: ysollw, ysolsw, ysollwdown REAL, DIMENSION(klon) :: yfder REAL, DIMENSION(klon) :: yrads,yrugoro REAL, DIMENSION(klon) :: yfluxlat REAL, DIMENSION(klon) :: y_d_ts REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1 REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q REAL, DIMENSION(klon) :: u1lay, v1lay REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m REAL, DIMENSION(klon) :: yustar REAL, DIMENSION(klon) :: yu10mx REAL, DIMENSION(klon) :: yu10my REAL, DIMENSION(klon) :: ywindsp REAL, DIMENSION(klon) :: yt10m, yq10m REAL, DIMENSION(klon) :: ypblh REAL, DIMENSION(klon) :: ylcl REAL, DIMENSION(klon) :: ycapCL REAL, DIMENSION(klon) :: yoliqCL REAL, DIMENSION(klon) :: ycteiCL REAL, DIMENSION(klon) :: ypblT REAL, DIMENSION(klon) :: ytherm REAL, DIMENSION(klon) :: ytrmb1 REAL, DIMENSION(klon) :: ytrmb2 REAL, DIMENSION(klon) :: ytrmb3 REAL, DIMENSION(klon) :: uzon, vmer REAL, DIMENSION(klon) :: tair1, qair1, tairsol REAL, DIMENSION(klon) :: psfce, patm REAL, DIMENSION(klon) :: qairsol, zgeo1 REAL, DIMENSION(klon) :: rugo1 REAL, DIMENSION(klon) :: yfluxsens, swdown REAL, DIMENSION(klon) :: petAcoef, peqAcoef, petBcoef, peqBcoef REAL, DIMENSION(klon) :: ypsref, epot_air REAL, DIMENSION(klon) :: yevap, ytsurf_new, yalb_new REAL, DIMENSION(klon) :: pctsrf_nsrf REAL, DIMENSION(klon) :: ztsol REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm REAL, DIMENSION(klon,klev) :: yu, yv REAL, DIMENSION(klon,klev) :: yt, yq REAL, DIMENSION(klon,klev) :: ypplay, ydelp REAL, DIMENSION(klon,klev) :: delp REAL, DIMENSION(klon,klev+1) :: ypaprs REAL, DIMENSION(klon,nsoilmx) :: ytsoil REAL, DIMENSION(klon,nbsrf) :: pctsrf_pot CHARACTER(len=80) :: abort_message CHARACTER(len=20) :: modname = 'pbl_surface' LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples LOGICAL, PARAMETER :: check=.FALSE. ! For debugging with IOIPSL INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg REAL :: zjulian REAL, DIMENSION(klon) :: tabindx REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat REAL, DIMENSION(iim,jjm+1) :: debugtab REAL, DIMENSION(klon,nbsrf) :: pblh REAL, DIMENSION(klon,nbsrf) :: plcl REAL, DIMENSION(klon,nbsrf) :: capCL REAL, DIMENSION(klon,nbsrf) :: oliqCL REAL, DIMENSION(klon,nbsrf) :: cteiCL REAL, DIMENSION(klon,nbsrf) :: pblT REAL, DIMENSION(klon,nbsrf) :: therm REAL, DIMENSION(klon,nbsrf) :: trmb1 REAL, DIMENSION(klon,nbsrf) :: trmb2 REAL, DIMENSION(klon,nbsrf) :: trmb3 REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m REAL, DIMENSION(klon,nbsrf) :: zx_qs1, zx_t1 REAL, DIMENSION(klon,nbsrf) :: zdelta1, zcor1 !jg+ temporary test REAL, DIMENSION(klon,klev) :: y_flux_u_old, y_flux_v_old REAL, DIMENSION(klon,klev) :: y_d_u_old, y_d_v_old !jg- !**************************************************************************************** ! End of declarations !**************************************************************************************** !**************************************************************************************** ! 1) Initialisation and validation tests ! Only done first time entering this subroutine ! !**************************************************************************************** IF (first_call) THEN first_call=.FALSE. ! Initilize debug IO IF (debugindex .AND. mpi_size==1) THEN ! initialize IOIPSL output idayref = day_ini CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian) CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon) DO i = 1, iim zx_lon(i,1) = rlon(i+1) zx_lon(i,jjm+1) = rlon(i+1) ENDDO CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat) CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), & 1,iim,1,jjm+1, & itau_phy,zjulian,dtime,nhoridbg,nidbg) ! no vertical axis cl_surf(1)='ter' cl_surf(2)='lic' cl_surf(3)='oce' cl_surf(4)='sic' DO nsrf=1,nbsrf CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, & jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime) END DO CALL histend(nidbg) CALL histsync(nidbg) END IF ENDIF !**************************************************************************************** ! 2) Initialization to zero ! Done for all local variables that will be compressed later ! and argument with INTENT(OUT) !**************************************************************************************** cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0 ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0 ; zu1 = 0.0 zv1 = 0.0 ; yqsurf = 0.0 ; yalb = 0.0 ; yalblw = 0.0 yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0 ysollw = 0.0 ; ysollwdown = 0.0 ; yrugos = 0.0 ; yu1 = 0.0 yv1 = 0.0 ; yrads = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0 ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0 yq = 0.0 ; pctsrf_new = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0 yrugoro = 0.0 ; yu10mx = 0.0 ; yu10my = 0.0 ; ywindsp = 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_u = 0.0 ; d_v = 0.0 ; zcoefh = 0.0 ; yqsol = 0.0 ytherm = 0.0 ytsoil = 999999. !**************************************************************************************** ! 3) - Calculate pressure thickness of each layer ! - Calculate the wind at first layer ! !**************************************************************************************** DO k = 1, klev DO i = 1, klon delp(i,k) = paprs(i,k)-paprs(i,k+1) ENDDO ENDDO DO i = 1, klon zx_alf1 = 1.0 zx_alf2 = 1.0 - zx_alf1 u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2 v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2 ENDDO !**************************************************************************************** ! Test for rugos........ from physiq.. A la fin plutot??? ! Calcul de l'abedo moyen par maille !**************************************************************************************** zxrugs(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015) zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO ! Calcul de l'abedo moyen par maille albsol(:) = 0.0 albsollw(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon albsol(i) = albsol(i) + albe(i,nsrf) * pctsrf(i,nsrf) albsollw(i) = albsollw(i) + alblw(i,nsrf) * pctsrf(i,nsrf) ENDDO ENDDO ! Calcule de ztsol (aussi fait dans physiq.F, pourrait etre un argument) ztsol(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO ! Repartition du longwave par sous-surface linearisee DO nsrf = 1, nbsrf DO i = 1, klon sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf)) solsw(i,nsrf) = solsw_m(i)*(1.-albe(i,nsrf))/(1.-albsol(i)) ENDDO ENDDO DO i = 1, klon sollwdown(i) = sollw_m(i) + RSIGMA*ztsol(i)**4 ENDDO !**************************************************************************************** ! 4) Loop over different surfaces ! ! All points with a possibility of the current surface are used. This is done ! to allow the sea-ice to appear or disappear. It is considered here that the ! entier domaine of the ocean possibly can contain sea-ice. ! !**************************************************************************************** pctsrf_pot = pctsrf pctsrf_pot(:,is_oce) = 1. - zmasq(:) pctsrf_pot(:,is_sic) = 1. - zmasq(:) loop_nbsrf: DO nsrf = 1, nbsrf ! Search for index(ni) and size(knon) of domaine to treat ni(:) = 0 knon = 0 DO i = 1, klon IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN knon = knon + 1 ni(knon) = i ENDIF ENDDO ! write index, with IOIPSL IF (debugindex .AND. mpi_size==1) THEN tabindx(:)=0. DO i=1,knon tabindx(i)=FLOAT(i) END DO debugtab(:,:) = 0. ndexbg(:) = 0 CALL gath2cpl(tabindx,debugtab,knon,ni) CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg) ENDIF !**************************************************************************************** ! 5) Compress variables ! !**************************************************************************************** DO j = 1, knon i = ni(j) ypct(j) = pctsrf(i,nsrf) yts(j) = ts(i,nsrf) ysnow(j) = snow(i,nsrf) yqsurf(j) = qsurf(i,nsrf) yalb(j) = albe(i,nsrf) yalblw(j) = alblw(i,nsrf) yrain_f(j) = rain_f(i) ysnow_f(j) = snow_f(i) yagesno(j) = agesno(i,nsrf) yfder(j) = fder(i) ysolsw(j) = solsw(i,nsrf) ysollw(j) = sollw(i,nsrf) ysollwdown(j) = sollwdown(i) yrugos(j) = rugos(i,nsrf) yrugoro(j) = rugoro(i) yu1(j) = u1lay(i) yv1(j) = v1lay(i) yrads(j) = ysolsw(j)+ ysollw(j) ypaprs(j,klev+1) = paprs(i,klev+1) yu10mx(j) = u10m(i,nsrf) yu10my(j) = v10m(i,nsrf) ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) ) END DO DO k = 1, klev DO j = 1, knon i = ni(j) ypaprs(j,k) = paprs(i,k) ypplay(j,k) = pplay(i,k) ydelp(j,k) = delp(i,k) yu(j,k) = u(i,k) yv(j,k) = v(i,k) yt(j,k) = t(i,k) yq(j,k) = q(i,k) ENDDO ENDDO DO k = 1, nsoilmx DO j = 1, knon i = ni(j) ytsoil(j,k) = ftsoil(i,k,nsrf) END DO END DO ! qsol(water height in soil) only for bucket continental model IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN DO j = 1, knon i = ni(j) yqsol(j) = qsol(i) END DO ENDIF !**************************************************************************************** ! 6) Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the ! atmosphere and coefficients for turbulent diffusion at surface(Cdrag). ! !**************************************************************************************** CALL coef_diff_turb(dtime, nsrf, knon, ni, & ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, & ycoefm, ycoefh) !jg+ !**************************************************************************************** ! => Old method ! Calculer la diffusion des vitesses "u" et "v" ! Output can be used : y_d_u_old, y_flux_u_old, y_d_v_old, y_flux_v_old ! !**************************************************************************************** CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp, & y_d_u_old, y_flux_u_old) CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp, & y_d_v_old, y_flux_v_old) ! save_flux est utile pour pouvoir utilise calcul_flux_vent plus tard CALL save_flux(klon, y_flux_u_old(:,1), y_flux_v_old(:,1)) !jg- !**************************************************************************************** ! ! 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. ! !**************************************************************************************** ! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, & ydelp, yt, yq, dtime, & petAcoef, peqAcoef, petBcoef, peqBcoef) ! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv) !**************************************************************************************** ! 9) Small calculations ! !**************************************************************************************** ypsref(:) = ypaprs(:,1) epot_air(:) = 0.0 epot_air(1:knon) = RCPD*yt(1:knon,1)*(ypsref(1:knon)/ypplay(1:knon,1))**RKAPPA swdown(:) = 0.0 IF (nsrf .EQ. is_ter) THEN swdown(1:knon) = ysolsw(1:knon)/(1-yalb(1:knon)) ELSE swdown(1:knon) = ysolsw(1:knon) ENDIF ! constant CO2 r_co2_ppm(:) = co2_ppm !**************************************************************************************** ! ! 10) Switch selon current surface ! It is necessary to start with the continental surfaces because the ocean ! needs their run-off. ! !**************************************************************************************** SELECT CASE(nsrf) CASE(is_ter) CALL surf_land(itap, dtime, date0, jour, knon, ni,& rlon, rlat, & debut, lafin, ydelp(:,1), epot_air, r_co2_ppm, ysollwdown, ysolsw, swdown, & yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),& petAcoef, peqAcoef, petBcoef, peqBcoef, & ypsref, yu1, yv1, yrugoro, pctsrf, & yrads, ysnow, yqsurf, yqsol, yagesno, & ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, & ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf) CASE(is_lic) CALL surf_landice(itap, dtime, knon, ni, & yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),& petAcoef, peqAcoef, petBcoef, peqBcoef, & ypsref, yu1, yv1, yrugoro, pctsrf, & yrads, ysnow, yqsurf, yqsol, yagesno, & ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, & ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf) CASE(is_oce) CALL surf_ocean(rlon, rlat, ysollw, yalb, & yrugos, ywindsp, rmu0, & yfder, & itap, dtime, jour, knon, ni, & debut, swdown, & ypplay(:,1), ycoefh(:,1), ycoefm(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),& petAcoef, peqAcoef, petBcoef, peqBcoef, & ypsref, yu1, yv1, yrugoro, pctsrf, & yrads, ysnow, yqsurf, yagesno, & yz0_new, yalblw, yevap, yfluxsens, yfluxlat, & ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf) CASE(is_sic) CALL surf_seaice( & rlon, rlat, ysollw, yalb, & yfder, & itap, dtime, jour, knon, ni, & debut, lafin, swdown, & yts, ypplay(:,1), ycoefh(:,1), yrain_f, ysnow_f, yt(:,1), yq(:,1),& petAcoef, peqAcoef, petBcoef, peqBcoef, & ypsref, yu1, yv1, yrugoro, pctsrf, & yrads, ysnow, yqsurf, yqsol, yagesno, & ytsoil, yz0_new, yalblw, yevap, yfluxsens, yfluxlat, & ytsurf_new, yalb_new, y_dflux_t, y_dflux_q, pctsrf_nsrf) CASE DEFAULT WRITE(lunout,*) 'Surface index = ', nsrf abort_message = 'Surface index not valid' CALL abort_gcm(modname,abort_message,1) END SELECT !**************************************************************************************** ! Save the fraction of this sub-surface ! !**************************************************************************************** pctsrf_new(:,nsrf) = pctsrf_nsrf(:) !**************************************************************************************** ! 11) - Calcul the increment of surface temperature ! - Update albedo ! !**************************************************************************************** y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon) yalb(1:knon) = yalb_new(1:knon) !**************************************************************************************** ! ! 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. ! !**************************************************************************************** ! H and Q y_flux_t1(:) = yfluxsens(:) y_flux_q1(:) = -yevap(:) CALL climb_hq_up(knon, dtime, yt, yq, & y_flux_q1, y_flux_t1, ypaprs, ypplay, & y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:)) ! U and V CALL climb_wind_up(knon, dtime, yu, yv, & y_flux_u, y_flux_v, y_d_u, y_d_v) !jg+ temporary for testing ! Use the results from old method IF (.NOT. newwind) THEN y_flux_u(:,:) = y_flux_u_old(:,:) y_flux_v(:,:) = y_flux_v_old(:,:) y_d_u(:,:) = y_d_u_old(:,:) y_d_v(:,:) = y_d_v_old(:,:) ENDIF !jg- DO j = 1, knon y_dflux_t(j) = y_dflux_t(j) * ypct(j) y_dflux_q(j) = y_dflux_q(j) * ypct(j) yu1(j) = yu1(j) * ypct(j) yv1(j) = yv1(j) * ypct(j) ENDDO !**************************************************************************************** ! 13) Transform variables for output format : ! - Decompress ! - Multiply with pourcentage of current surface ! - Cumulate in global variable ! !**************************************************************************************** DO k = 1, klev DO j = 1, knon i = ni(j) ycoefh(j,k) = ycoefh(j,k) * ypct(j) ycoefm(j,k) = ycoefm(j,k) * ypct(j) y_d_t(j,k) = y_d_t(j,k) * ypct(j) y_d_q(j,k) = y_d_q(j,k) * ypct(j) y_d_u(j,k) = y_d_u(j,k) * ypct(j) y_d_v(j,k) = y_d_v(j,k) * ypct(j) flux_t(i,k,nsrf) = y_flux_t(j,k) flux_q(i,k,nsrf) = y_flux_q(j,k) flux_u(i,k,nsrf) = y_flux_u(j,k) flux_v(i,k,nsrf) = y_flux_v(j,k) ENDDO ENDDO evap(:,nsrf) = - flux_q(:,1,nsrf) albe(:, nsrf) = 0. alblw(:, nsrf) = 0. snow(:, nsrf) = 0. qsurf(:, nsrf) = 0. rugos(:, nsrf) = 0. fluxlat(:,nsrf) = 0. DO j = 1, knon i = ni(j) d_ts(i,nsrf) = y_d_ts(j) albe(i,nsrf) = yalb(j) alblw(i,nsrf) = yalblw(j) snow(i,nsrf) = ysnow(j) qsurf(i,nsrf) = yqsurf(j) rugos(i,nsrf) = yz0_new(j) fluxlat(i,nsrf) = yfluxlat(j) agesno(i,nsrf) = yagesno(j) cdragh(i) = cdragh(i) + ycoefh(j,1) cdragm(i) = cdragm(i) + ycoefm(j,1) dflux_t(i) = dflux_t(i) + y_dflux_t(j) dflux_q(i) = dflux_q(i) + y_dflux_q(j) zu1(i) = zu1(i) + yu1(j) zv1(i) = zv1(i) + yv1(j) END DO IF ( nsrf .EQ. is_ter ) THEN DO j = 1, knon i = ni(j) qsol(i) = yqsol(j) END DO END IF ftsoil(:,:,nsrf) = 0. DO k = 1, nsoilmx DO j = 1, knon i = ni(j) ftsoil(i, k, nsrf) = ytsoil(j,k) END DO END DO #ifdef CRAY DO k = 1, klev DO j = 1, knon i = ni(j) #else DO j = 1, knon i = ni(j) DO k = 1, klev #endif d_t(i,k) = d_t(i,k) + y_d_t(j,k) d_q(i,k) = d_q(i,k) + y_d_q(j,k) d_u(i,k) = d_u(i,k) + y_d_u(j,k) d_v(i,k) = d_v(i,k) + y_d_v(j,k) zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k) #ifdef CRAY END DO END DO #else END DO END DO #endif !**************************************************************************************** ! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m ! Call HBTM ! !**************************************************************************************** t2m(:,nsrf) = 0. q2m(:,nsrf) = 0. u10m(:,nsrf) = 0. v10m(:,nsrf) = 0. pblh(:,nsrf) = 0. ! Hauteur de couche limite plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA capCL(:,nsrf) = 0. ! CAPE de couche limite oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite therm(:,nsrf) = 0. trmb1(:,nsrf) = 0. ! deep_cape trmb2(:,nsrf) = 0. ! inhibition trmb3(:,nsrf) = 0. ! Point Omega #undef T2m #define T2m #ifdef T2m ! diagnostic t,q a 2m et u, v a 10m DO j=1, knon i = ni(j) uzon(j) = yu(j,1) + y_d_u(j,1) vmer(j) = yv(j,1) + y_d_v(j,1) tair1(j) = yt(j,1) + y_d_t(j,1) qair1(j) = yq(j,1) + y_d_q(j,1) zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) & * (ypaprs(j,1)-ypplay(j,1)) tairsol(j) = yts(j) + y_d_ts(j) rugo1(j) = yrugos(j) IF(nsrf.EQ.is_oce) THEN rugo1(j) = rugos(i,nsrf) ENDIF psfce(j)=ypaprs(j,1) patm(j)=ypplay(j,1) qairsol(j) = yqsurf(j) END DO ! Calculate the temperature et relative humidity at 2m and the wind at 10m CALL stdlevvar(klon, knon, nsrf, zxli, & uzon, vmer, tair1, qair1, zgeo1, & tairsol, qairsol, rugo1, psfce, patm, & yt2m, yq2m, yt10m, yq10m, yu10m, yustar) DO j=1, knon i = ni(j) t2m(i,nsrf)=yt2m(j) q2m(i,nsrf)=yq2m(j) ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2) v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2) END DO CALL HBTM(knon, ypaprs, ypplay, & yt2m,yt10m,yq2m,yq10m,yustar, & y_flux_t,y_flux_q,yu,yv,yt,yq, & ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, & ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl) DO j=1, knon i = ni(j) pblh(i,nsrf) = ypblh(j) plcl(i,nsrf) = ylcl(j) capCL(i,nsrf) = ycapCL(j) oliqCL(i,nsrf) = yoliqCL(j) cteiCL(i,nsrf) = ycteiCL(j) pblT(i,nsrf) = ypblT(j) therm(i,nsrf) = ytherm(j) trmb1(i,nsrf) = ytrmb1(j) trmb2(i,nsrf) = ytrmb2(j) trmb3(i,nsrf) = ytrmb3(j) END DO #else ! not defined T2m ! No calculation ! Set output variables to zero DO j = 1, knon i = ni(j) pblh(i,nsrf) = 0. plcl(i,nsrf) = 0. capCL(i,nsrf) = 0. oliqCL(i,nsrf) = 0. cteiCL(i,nsrf) = 0. pblT(i,nsrf) = 0. therm(i,nsrf) = 0. trmb1(i,nsrf) = 0. trmb2(i,nsrf) = 0. trmb3(i,nsrf) = 0. END DO DO j = 1, knon i = ni(j) t2m(i,nsrf)=0. q2m(i,nsrf)=0. u10m(i,nsrf)=0. v10m(i,nsrf)=0. END DO #endif !**************************************************************************************** ! 15) End of loop over different surfaces ! !**************************************************************************************** END DO loop_nbsrf !**************************************************************************************** ! 16) Calculate the mean value over all sub-surfaces for som variables ! ! NB!!! jg : Pour garder la convergence numerique j'utilise pctsrf_new comme c'etait ! fait dans physiq.F mais ca devrait plutot etre pctsrf???!!!!! A verifier! !**************************************************************************************** zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0 zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0 DO nsrf = 1, nbsrf DO k = 1, klev DO i = 1, klon zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf_new(i,nsrf) zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf_new(i,nsrf) zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf_new(i,nsrf) zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf_new(i,nsrf) END DO END DO END DO DO i = 1, klon zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i) ENDDO DO i = 1, klon IF ( ABS( pctsrf_new(i, is_ter) + pctsrf_new(i, is_lic) + & pctsrf_new(i, is_oce) + pctsrf_new(i, is_sic) - 1.) .GT. EPSFRA) & THEN WRITE(*,*) 'physiq : pb sous surface au point ', i, & pctsrf_new(i, 1 : nbsrf) ENDIF ENDDO ! ! Incrementer la temperature du sol ! zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0 zt2m(:) = 0.0 ; zq2m(:) = 0.0 zu10m(:) = 0.0 ; zv10m(:) = 0.0 s_pblh(:) = 0.0 ; s_plcl(:) = 0.0 s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0 s_cteiCL(:) = 0.0; s_pblT(:) = 0.0 s_therm(:) = 0.0 ; s_trmb1(:) = 0.0 s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf) wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) & + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf_new(i,nsrf) wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * & pctsrf_new(i,nsrf) zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf_new(i,nsrf) zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf_new(i,nsrf) zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf_new(i,nsrf) zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf_new(i,nsrf) zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf_new(i,nsrf) zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf_new(i,nsrf) s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf_new(i,nsrf) s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf_new(i,nsrf) s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf_new(i,nsrf) s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf_new(i,nsrf) s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf_new(i,nsrf) s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf_new(i,nsrf) s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf_new(i,nsrf) s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf_new(i,nsrf) s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf_new(i,nsrf) s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf_new(i,nsrf) END DO END DO IF (check) THEN amn=MIN(ts(1,is_ter),1000.) amx=MAX(ts(1,is_ter),-1000.) DO i=2, klon amn=MIN(ts(i,is_ter),amn) amx=MAX(ts(i,is_ter),amx) ENDDO PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx ENDIF ! ! 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 DO i = 1, klon fder(i) = - 4.0*RSIGMA*zxtsol(i)**3 ENDDO zxqsurf(:) = 0.0 zxsnow(:) = 0.0 DO nsrf = 1, nbsrf DO i = 1, klon zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf_new(i,nsrf) zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf_new(i,nsrf) END DO END DO ! !IM Calculer l'humidite relative a 2m (rh2m) pour diagnostique !IM ajout dependance type surface rh2m(:) = 0.0 qsat2m(:) = 0.0 DO i = 1, klon DO nsrf=1, nbsrf zx_t1(i,nsrf) = t2m(i,nsrf) IF (thermcep) THEN zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf))) zx_qs1(i,nsrf) = r2es * & FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1) zx_qs1(i,nsrf) = MIN(0.5,zx_qs1(i,nsrf)) zcor1(i,nsrf) = 1./(1.-retv*zx_qs1(i,nsrf)) zx_qs1(i,nsrf) = zx_qs1(i,nsrf)*zcor1(i,nsrf) END IF zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf) zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf) rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf_new(i,nsrf) qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf_new(i,nsrf) END DO END DO ! Some of the module declared variables are returned for printing in physiq.F qsol_d(:) = qsol(:) evap_d(:,:) = evap(:,:) rugos_d(:,:) = rugos(:,:) agesno_d(:,:) = agesno(:,:) END SUBROUTINE pbl_surface ! !**************************************************************************************** ! SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, & evap_rst, rugos_rst, agesno_rst, ftsoil_rst) INCLUDE "indicesol.inc" INCLUDE "dimsoil.h" ! Ouput variables !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst !**************************************************************************************** ! Return module variables for writing to restart file ! !**************************************************************************************** qsol_rst(:) = qsol(:) fder_rst(:) = fder(:) snow_rst(:,:) = snow(:,:) qsurf_rst(:,:) = qsurf(:,:) evap_rst(:,:) = evap(:,:) rugos_rst(:,:) = rugos(:,:) agesno_rst(:,:) = agesno(:,:) ftsoil_rst(:,:,:) = ftsoil(:,:,:) !**************************************************************************************** ! Deallocate module variables ! !**************************************************************************************** DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil) END SUBROUTINE pbl_surface_final ! !**************************************************************************************** ! END MODULE pbl_surface_mod