c c $Header$ c SUBROUTINE physiq (nlon,nlev,nqmax, . debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys, . paprs,pplay,pphi,pphis,paire,presnivs,clesphy0, . u,v,t,qx, . omega, cufi, cvfi, . d_u, d_v, d_t, d_qx, d_ps) USE ioipsl USE histcom USE writephys IMPLICIT none c====================================================================== c c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 c c Objet: Moniteur general de la physique du modele cAA Modifications quant aux traceurs : cAA - uniformisation des parametrisations ds phytrac cAA - stockage des moyennes des champs necessaires cAA en mode traceur off-line c====================================================================== c modif ( P. Le Van , 12/10/98 ) c c Arguments: c c nlon----input-I-nombre de points horizontaux c nlev----input-I-nombre de couches verticales c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1 c debut---input-L-variable logique indiquant le premier passage c lafin---input-L-variable logique indiquant le dernier passage c rjour---input-R-numero du jour de l'experience c gmtime--input-R-temps universel dans la journee (0 a 86400 s) c pdtphys-input-R-pas d'integration pour la physique (seconde) c paprs---input-R-pression pour chaque inter-couche (en Pa) c pplay---input-R-pression pour le mileu de chaque couche (en Pa) c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol) c pphis---input-R-geopotentiel du sol c paire---input-R-aire de chaque maille c presnivs-input_R_pressions approximat. des milieux couches ( en PA) c u-------input-R-vitesse dans la direction X (de O a E) en m/s c v-------input-R-vitesse Y (de S a N) en m/s c t-------input-R-temperature (K) c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs c d_t_dyn-input-R-tendance dynamique pour "t" (K/s) c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s) c omega---input-R-vitesse verticale en Pa/s c cufi----input-R-resolution des mailles en x (m) c cvfi----input-R-resolution des mailles en y (m) c c d_u-----output-R-tendance physique de "u" (m/s/s) c d_v-----output-R-tendance physique de "v" (m/s/s) c d_t-----output-R-tendance physique de "t" (K/s) c d_qx----output-R-tendance physique de "qx" (kg/kg/s) c d_ps----output-R-tendance physique de la pression au sol c====================================================================== #include "dimensions.h" integer jjmp1 parameter (jjmp1=jjm+1-1/jjm) #include "dimphy.h" #include "regdim.h" #include "indicesol.h" #include "dimsoil.h" #include "clesphys.h" #include "control.h" #include "temps.h" c====================================================================== LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE PARAMETER (ok_cvl=.TRUE.) LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface PARAMETER (ok_gust=.FALSE.) c====================================================================== LOGICAL check ! Verifier la conservation du modele en eau PARAMETER (check=.FALSE.) LOGICAL ok_stratus ! Ajouter artificiellement les stratus PARAMETER (ok_stratus=.FALSE.) c====================================================================== c Parametres lies au coupleur OASIS: #include "oasis.h" INTEGER,SAVE :: npas, nexca logical rnpb parameter(rnpb=.true.) c PARAMETER (npas=1440) c PARAMETER (nexca=48) EXTERNAL fromcpl, intocpl, inicma c ocean = type de modele ocean a utiliser: force, slab, couple character*6 ocean SAVE ocean c parameter (ocean = 'force ') c parameter (ocean = 'couple') logical ok_ocean c====================================================================== c Clef controlant l'activation du cycle diurne: ccc LOGICAL cycle_diurne ccc PARAMETER (cycle_diurne=.FALSE.) c====================================================================== c Modele thermique du sol, a activer pour le cycle diurne: ccc LOGICAL soil_model ccc PARAMETER (soil_model=.FALSE.) logical ok_veget save ok_veget c parameter (ok_veget = .true.) c parameter (ok_veget = .false.) c====================================================================== c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans c le calcul du rayonnement est celle apres la precipitation des nuages. c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre c la condensation et la precipitation. Cette cle augmente les impacts c radiatifs des nuages. ccc LOGICAL new_oliq ccc PARAMETER (new_oliq=.FALSE.) c====================================================================== c Clefs controlant deux parametrisations de l'orographie: cc LOGICAL ok_orodr ccc PARAMETER (ok_orodr=.FALSE.) ccc LOGICAL ok_orolf ccc PARAMETER (ok_orolf=.FALSE.) c====================================================================== LOGICAL ok_journe ! sortir le fichier journalier save ok_journe c PARAMETER (ok_journe=.true.) integer lev_histday save lev_histday data lev_histday/1/ c LOGICAL ok_mensuel ! sortir le fichier mensuel save ok_mensuel c PARAMETER (ok_mensuel=.true.) c LOGICAL ok_mensuelNMC ! sortir le fichier mensuel niveaux NMC PARAMETER (ok_mensuelNMC=.true.) c save ok_mensuelNMC c LOGICAL ok_instan ! sortir le fichier instantane save ok_instan c PARAMETER (ok_instan=.true.) c LOGICAL ok_region ! sortir le fichier regional PARAMETER (ok_region=.FALSE.) c====================================================================== c INTEGER ivap ! indice de traceurs pour vapeur d'eau PARAMETER (ivap=1) INTEGER iliq ! indice de traceurs pour eau liquide PARAMETER (iliq=2) c c c Variables argument: c INTEGER nlon INTEGER nlev INTEGER nqmax REAL rjourvrai, rjour_ecri REAL gmtime REAL pdtphys LOGICAL debut, lafin REAL paprs(klon,klev+1) REAL pplay(klon,klev) REAL pphi(klon,klev) REAL pphis(klon) REAL paire(klon) REAL presnivs(klev) REAL znivsig(klev) REAL zsurf(nbsrf) real cufi(klon), cvfi(klon) REAL u(klon,klev) REAL v(klon,klev) REAL t(klon,klev) REAL qx(klon,klev,nqmax) REAL t_ancien(klon,klev), q_ancien(klon,klev) SAVE t_ancien, q_ancien LOGICAL ancien_ok SAVE ancien_ok REAL d_t_dyn(klon,klev) REAL d_q_dyn(klon,klev) REAL omega(klon,klev) REAL d_u(klon,klev) REAL d_v(klon,klev) REAL d_t(klon,klev) REAL d_qx(klon,klev,nqmax) REAL d_ps(klon) INTEGER klevp1, klevm1 PARAMETER(klevp1=klev+1,klevm1=klev-1) #include "raddim.h" c REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2) SAVE swdn0 , swdn, swup0, swup c vents meridien et zonal a un niveau de pression real u1000(klon), v1000(klon) !vents a 1000 hPa real u925(klon), v925(klon) !vents a 925 hPa real u850(klon),v850(klon) !vents a 850 hPa real u700(klon),v700(klon) real u600(klon),v600(klon) real u500(klon),v500(klon) real u400(klon),v400(klon) real u300(klon),v300(klon) real u250(klon),v250(klon) real u200(klon),v200(klon) real u150(klon),v150(klon) real u100(klon),v100(klon) real u70(klon),v70(klon) real u50(klon),v50(klon) real u30(klon),v30(klon) real u20(klon),v20(klon) real u10(klon),v10(klon) real phi500(klon),w500(klon) c prw: precipitable water real prw(klon) REAL convliq(klon,klev) ! eau liquide nuageuse convective REAL convfra(klon,klev) ! fraction nuageuse convective REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree INTEGER linv, kp1 c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2) c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg) REAL flwp(klon), fiwp(klon) REAL flwc(klon,klev), fiwc(klon,klev) REAL flwp_c(klon), fiwp_c(klon) REAL flwc_c(klon,klev), fiwc_c(klon,klev) REAL flwp_s(klon), fiwp_s(klon) REAL flwc_s(klon,klev), fiwc_s(klon,klev) c ISCCP simulator v3.4 c dans clesphys.h top_height, overlap cv3.4 INTEGER debug, debugcol INTEGER npoints PARAMETER(npoints=klon) c INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night INTEGER nregISCtot PARAMETER(nregISCtot=1) c c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire c y compris pour 1 point c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude) c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude) INTEGER imin_debut, nbpti INTEGER jmin_debut, nbptj c REAL nbsunlit(nregISCtot,klon) !nbsunlit : moyenne de sunlit INTEGER ncol, seed(klon) c ncol = nb. de sous-colonnes pour chaque maille du GCM c PARAMETER(ncol=100) c PARAMETER(ncol=625) c PARAMETER(ncol=10) PARAMETER(ncol=25) REAL tautab(0:255) INTEGER invtau(-20:45000) REAL emsfc_lw PARAMETER(emsfc_lw=0.99) REAL ran0 ! type for random number fuction c REAL cldtot(klon,klev) c variables de haut en bas pour le simulateur ISCCP REAL dtau_s(klon,klev) !tau nuages startiformes REAL dtau_c(klon,klev) !tau nuages convectifs REAL dem_s(klon,klev) !emissivite nuages startiformes REAL dem_c(klon,klev) !emissivite nuages convectifs c c variables de haut en bas pour le simulateur ISCCP REAL pfull(klon,klev) REAL phalf(klon,klev+1) REAL qv(klon,klev) REAL cc(klon,klev) REAL conv(klon,klev) REAL dtau_sH2B(klon,klev) REAL dtau_cH2B(klon,klev) REAL at(klon,klev) REAL dem_sH2B(klon,klev) REAL dem_cH2B(klon,klev) c output from ISCCP simulator REAL fq_isccp(klon,7,7) REAL totalcldarea(klon) REAL meanptop(klon) REAL meantaucld(klon) REAL boxtau(klon,ncol) REAL boxptop(klon,ncol) c INTEGER l, ni, nj, kmax, lmax PARAMETER(kmax=8, lmax=8) INTEGER kmaxm1, lmaxm1 PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1) INTEGER iimx7, jjmx7, jjmp1x7 PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1, .jjmp1x7=jjmp1*lmaxm1) REAL fq4d(iim,jjmp1,kmaxm1,lmaxm1) REAL fq3d(iimx7, jjmp1x7) c INTEGER iw, iwmax REAL wmin, pas_w c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30) PARAMETER(wmin=-200.,pas_w=10.,iwmax=40) REAL o500(klon) c cIM: nbregdyn = nbre regions pour calculs statistiques sur output du ISCCP cIM: dynamiques INTEGER nreg, nbregdyn PARAMETER(nbregdyn=5) REAL histoW(kmaxm1,lmaxm1,iwmax,nbregdyn) REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbregdyn) REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbregdyn) SAVE nhistoWt INTEGER pct_ocean(klon,nbregdyn) REAL rlonPOS(klon) c sorties ISCCP logical ok_isccp real ecrit_isccp integer nid_isccp save ok_isccp, ecrit_isccp, nid_isccp #define histISCCP #undef histISCCP #ifdef histISCCP c data ok_isccp,ecrit_isccp/.true.,0.125/ c data ok_isccp,ecrit_isccp/.true.,1./ data ok_isccp/.true./ #else data ok_isccp/.false./ #endif c sorties statistiques regime dynamique logical ok_regdyn real ecrit_regdyn integer nid_regdyn save ok_regdyn, ecrit_regdyn, nid_regdyn #undef histREGDYN #define histREGDYN #ifdef histREGDYN c data ok_regdyn,ecrit_regdyn/.true.,0.125/ c data ok_regdyn,ecrit_regdyn/.true.,1./ data ok_regdyn/.true./ #else data ok_regdyn/.false./ #endif REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax) DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./ DATA zx_pc/50., 180., 310., 440., 560., 680., 800./ c cldtopres pression au sommet des nuages REAL cldtopres(lmaxm1) DATA cldtopres/50., 180., 310., 440., 560., 680., 800./ INTEGER komega, nhoriRD c taulev: numero du niveau de tau dans les sorties ISCCP CHARACTER *4 taulev(kmaxm1) DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/ REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7) INTEGER nhorix7 cIM: region='3d' <==> sorties en global CHARACTER*3 region PARAMETER(region='3d') c logical ok_hf real ecrit_hf integer nid_hf save ok_hf, ecrit_hf, nid_hf c QUESTION : noms de variables ? #define histhf #ifdef histhf data ok_hf,ecrit_hf/.true.,0.25/ #else data ok_hf/.false./ #endif INTEGER longcles PARAMETER ( longcles = 20 ) REAL clesphy0( longcles ) c c Variables quasi-arguments c REAL xjour SAVE xjour c c c Variables propres a la physique c REAL dtime SAVE dtime ! pas temporel de la physique c INTEGER radpas SAVE radpas ! frequence d'appel rayonnement c REAL radsol(klon) SAVE radsol ! bilan radiatif au sol calcule par code radiatif c REAL rlat(klon) SAVE rlat ! latitude pour chaque point c REAL rlon(klon) SAVE rlon ! longitude pour chaque point c cc INTEGER iflag_con cc SAVE iflag_con ! indicateur de la convection c INTEGER itap SAVE itap ! compteur pour la physique c REAL co2_ppm_etat0 c REAL solaire_etat0 c real slp(klon) ! sea level pressure REAL ftsol(klon,nbsrf) SAVE ftsol ! temperature du sol c REAL ftsoil(klon,nsoilmx,nbsrf) SAVE ftsoil ! temperature dans le sol c REAL fevap(klon,nbsrf) SAVE fevap ! evaporation REAL fluxlat(klon,nbsrf) SAVE fluxlat c REAL deltat(klon) SAVE deltat ! ecart avec la SST de reference c REAL fqsurf(klon,nbsrf) SAVE fqsurf ! humidite de l'air au contact de la surface c REAL qsol(klon) SAVE qsol ! hauteur d'eau dans le sol c REAL fsnow(klon,nbsrf) SAVE fsnow ! epaisseur neigeuse c REAL falbe(klon,nbsrf) SAVE falbe ! albedo par type de surface REAL falblw(klon,nbsrf) SAVE falblw ! albedo par type de surface c c c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM): c REAL zmea(klon) SAVE zmea ! orographie moyenne c REAL zstd(klon) SAVE zstd ! deviation standard de l'OESM c REAL zsig(klon) SAVE zsig ! pente de l'OESM c REAL zgam(klon) save zgam ! anisotropie de l'OESM c REAL zthe(klon) SAVE zthe ! orientation de l'OESM c REAL zpic(klon) SAVE zpic ! Maximum de l'OESM c REAL zval(klon) SAVE zval ! Minimum de l'OESM c REAL rugoro(klon) SAVE rugoro ! longueur de rugosite de l'OESM c REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon) c REAL zuthe(klon),zvthe(klon) SAVE zuthe SAVE zvthe INTEGER igwd,idx(klon),itest(klon) c REAL agesno(klon,nbsrf) SAVE agesno ! age de la neige c REAL alb_neig(klon) SAVE alb_neig ! albedo de la neige cKE43 c Variables liees a la convection de K. Emanuel (sb): c REAL ema_workcbmf(klon) ! cloud base mass flux SAVE ema_workcbmf REAL ema_cbmf(klon) ! cloud base mass flux SAVE ema_cbmf REAL ema_pcb(klon) ! cloud base pressure SAVE ema_pcb REAL ema_pct(klon) ! cloud top pressure SAVE ema_pct REAL bas, top ! cloud base and top levels SAVE bas SAVE top REAL Ma(klon,klev) ! undilute upward mass flux SAVE Ma REAL qcondc(klon,klev) ! in-cld water content from convect SAVE qcondc REAL ema_work1(klon, klev), ema_work2(klon, klev) SAVE ema_work1, ema_work2 REAL wdn(klon), tdn(klon), qdn(klon) REAL wd(klon) ! sb SAVE wd ! sb c Variables locales pour la couche limite (al1): c cAl1 REAL pblh(klon) ! Hauteur de couche limite cAl1 SAVE pblh c34EK c c Variables locales: c REAL cdragh(klon) ! drag coefficient pour T and Q REAL cdragm(klon) ! drag coefficient pour vent cAA cAA Pour phytrac cAA REAL ycoefh(klon,klev) ! coef d'echange pour phytrac REAL yu1(klon) ! vents dans la premiere couche U REAL yv1(klon) ! vents dans la premiere couche V REAL ffonte(klon,nbsrf) !Flux thermique utilise pour fondre la neige REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface c !et necessaire pour limiter la c !hauteur de neige, en kg/m2/s REAL zxffonte(klon), zxfqcalving(klon) LOGICAL offline ! Controle du stockage ds "physique" PARAMETER (offline=.false.) INTEGER physid REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction save pfrac_impa REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation save pfrac_nucl REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1) save pfrac_1nucl REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction) REAL frac_nucl(klon,klev) ! idem (nucleation) cAA REAL rain_fall(klon) ! pluie REAL snow_fall(klon) ! neige save snow_fall, rain_fall REAL evap(klon), devap(klon) ! evaporation et sa derivee REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee REAL dlw(klon) ! derivee infra rouge REAL bils(klon) ! bilan de chaleur au sol REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque C type de sous-surface et pondere par la fraction REAL fder(klon) ! Derive de flux (sensible et latente) save fder REAL ve(klon) ! integr. verticale du transport meri. de l'energie REAL vq(klon) ! integr. verticale du transport meri. de l'eau REAL ue(klon) ! integr. verticale du transport zonal de l'energie REAL uq(klon) ! integr. verticale du transport zonal de l'eau c REAL frugs(klon,nbsrf) ! longueur de rugosite save frugs REAL zxrugs(klon) ! longueur de rugosite c c Conditions aux limites c INTEGER julien c INTEGER lmt_pas SAVE lmt_pas ! frequence de mise a jour REAL pctsrf(klon,nbsrf) SAVE pctsrf ! sous-fraction du sol REAL albsol(klon) SAVE albsol ! albedo du sol total REAL albsollw(klon) SAVE albsollw ! albedo du sol total REAL wo(klon,klev) SAVE wo ! ozone c====================================================================== c c Declaration des procedures appelees c EXTERNAL angle ! calculer angle zenithal du soleil EXTERNAL alboc ! calculer l'albedo sur ocean EXTERNAL ajsec ! ajustement sec EXTERNAL clmain ! couche limite EXTERNAL condsurf ! lire les conditions aux limites EXTERNAL conlmd ! convection (schema LMD) cKE43 EXTERNAL conema3 ! convect4.3 EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie) cAA EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie) c ! stockage des coefficients necessaires au c ! lessivage OFF-LINE et ON-LINE cAA EXTERNAL hgardfou ! verifier les temperatures EXTERNAL nuage ! calculer les proprietes radiatives EXTERNAL o3cm ! initialiser l'ozone EXTERNAL orbite ! calculer l'orbite terrestre EXTERNAL ozonecm ! prescrire l'ozone EXTERNAL phyetat0 ! lire l'etat initial de la physique EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique EXTERNAL radlwsw ! rayonnements solaire et infrarouge EXTERNAL suphec ! initialiser certaines constantes EXTERNAL transp ! transport total de l'eau et de l'energie EXTERNAL ecribina ! ecrire le fichier binaire global EXTERNAL ecribins ! ecrire le fichier binaire global EXTERNAL ecrirega ! ecrire le fichier binaire regional EXTERNAL ecriregs ! ecrire le fichier binaire regional cIM EXTERNAL haut2bas !variables de haut en bas c c Variables locales c real clwcon(klon,klev),rnebcon(klon,klev) real clwcon0(klon,klev),rnebcon0(klon,klev) save rnebcon, clwcon REAL rhcl(klon,klev) ! humiditi relative ciel clair REAL dialiq(klon,klev) ! eau liquide nuageuse REAL diafra(klon,klev) ! fraction nuageuse REAL cldliq(klon,klev) ! eau liquide nuageuse REAL cldfra(klon,klev) ! fraction nuageuse REAL cldtau(klon,klev) ! epaisseur optique REAL cldemi(klon,klev) ! emissivite infrarouge c CXXX PB REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v c REAL zxfluxt(klon, klev) REAL zxfluxq(klon, klev) REAL zxfluxu(klon, klev) REAL zxfluxv(klon, klev) CXXX REAL heat(klon,klev) ! chauffage solaire REAL heat0(klon,klev) ! chauffage solaire ciel clair REAL cool(klon,klev) ! refroidissement infrarouge REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon) real sollwdown(klon) ! downward LW flux at surface REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon) REAL albpla(klon) cIM cf. JLD REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface c Le rayonnement n'est pas calcule tous les pas, il faut donc c sauvegarder les sorties du rayonnement SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0 cccIM INTEGER itaprad SAVE itaprad c REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s) REAL conv_t(klon,klev) ! convergence de la temperature(K/s) c REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree c REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon) c REAL dist, rmu0(klon), fract(klon) REAL zdtime, zlongi c CHARACTER*2 str2 CHARACTER*2 iqn c REAL qcheck REAL z_avant(klon), z_apres(klon), z_factor(klon) LOGICAL zx_ajustq c REAL za, zb REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp real zqsat(klon,klev) INTEGER i, k, iq, ig, j, nsrf, ll REAL t_coup PARAMETER (t_coup=234.0) c REAL zphi(klon,klev) REAL zx_tmp_x(iim), zx_tmp_yjjmp1 REAL zx_relief(iim,jjmp1) REAL zx_aire(iim,jjmp1) cKE43 c Variables locales pour la convection de K. Emanuel (sb): c REAL upwd(klon,klev) ! saturated updraft mass flux REAL dnwd(klon,klev) ! saturated downdraft mass flux REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux REAL tvp(klon,klev) ! virtual temp of lifted parcel REAL cape(klon) ! CAPE SAVE cape CHARACTER*40 capemaxcels !max(CAPE) REAL pbase(klon) ! cloud base pressure SAVE pbase REAL bbase(klon) ! cloud base buoyancy SAVE bbase REAL rflag(klon) ! flag fonctionnement de convect INTEGER iflagctrl(klon) ! flag fonctionnement de convect c -- convect43: INTEGER ntra ! nb traceurs pour convect4.3 REAL pori_con(klon) ! pressure at the origin level of lifted parcel REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon) REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev) REAL dplcldt(klon), dplcldr(klon) c? . condm_con(klon,klev),conda_con(klon,klev), c? . mr_con(klon,klev),ep_con(klon,klev) c? . ,sadiab(klon,klev),wadiab(klon,klev) c -- c34EK c c Variables du changement c c con: convection c lsc: condensation a grande echelle (Large-Scale-Condensation) c ajs: ajustement sec c eva: evaporation de l'eau liquide nuageuse c vdf: couche limite (Vertical DiFfusion) REAL d_t_con(klon,klev),d_q_con(klon,klev) REAL d_u_con(klon,klev),d_v_con(klon,klev) REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev) REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev) REAL d_t_eva(klon,klev),d_q_eva(klon,klev) REAL rneb(klon,klev) c REAL pmfu(klon,klev), pmfd(klon,klev) REAL pen_u(klon,klev), pen_d(klon,klev) REAL pde_u(klon,klev), pde_d(klon,klev) INTEGER kcbot(klon), kctop(klon), kdtop(klon) REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) REAL prfl(klon,klev+1), psfl(klon,klev+1) c INTEGER ibas_con(klon), itop_con(klon) REAL rain_con(klon), rain_lsc(klon) REAL snow_con(klon), snow_lsc(klon) REAL d_ts(klon,nbsrf) c REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev) REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev) c REAL d_u_oro(klon,klev), d_v_oro(klon,klev) REAL d_t_oro(klon,klev) REAL d_u_lif(klon,klev), d_v_lif(klon,klev) REAL d_t_lif(klon,klev) REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev) real ratqsbas,ratqshaut save ratqsbas,ratqshaut, ratqs real zpt_conv(klon,klev) c Parametres lies au nouveau schema de nuages (SB, PDF) real fact_cldcon real facttemps logical ok_newmicro save ok_newmicro save fact_cldcon,facttemps real facteur integer iflag_cldcon save iflag_cldcon logical ptconv(klon,klev) c c Variables liees a l'ecriture de la bande histoire physique c INTEGER ecrit_mth SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel) c INTEGER ecrit_day SAVE ecrit_day ! frequence d'ecriture (fichier journalier) c INTEGER ecrit_ins SAVE ecrit_ins ! frequence d'ecriture (fichier instantane) c INTEGER ecrit_reg SAVE ecrit_reg ! frequence d'ecriture c integer itau_w ! pas de temps ecriture = itap + itau_phy c c c Variables locales pour effectuer les appels en serie c REAL t_seri(klon,klev), q_seri(klon,klev) REAL ql_seri(klon,klev),qs_seri(klon,klev) REAL u_seri(klon,klev), v_seri(klon,klev) c REAL tr_seri(klon,klev,nbtr) REAL d_tr(klon,klev,nbtr) REAL zx_rh(klon,klev) INTEGER length PARAMETER ( length = 100 ) REAL tabcntr0( length ) c INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev) REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev) REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1) c INTEGER nid_day, nid_mth, nid_ins, nid_nmc SAVE nid_day, nid_mth, nid_ins, nid_nmc c INTEGER nhori, nvert REAL zsto, zout real zjulian save zjulian character*20 modname character*80 abort_message logical ok_sync real date0 integer idayref C essai writephys integer fid_day, fid_mth, fid_ins parameter (fid_ins = 1, fid_day = 2, fid_mth = 3) integer prof2d_on, prof3d_on, prof2d_av, prof3d_av parameter (prof2d_on = 1, prof3d_on = 2, . prof2d_av = 3, prof3d_av = 4) character*30 nom_fichier character*10 varname character*40 vartitle character*20 varunits C Variables liees au bilan d'energie et d'enthalpi REAL ztsol(klon) REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot $ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec REAL d_h_vcol_phy REAL fs_bound, fq_bound SAVE d_h_vcol_phy REAL zero_v(klon) CHARACTER*15 ztit INTEGER ip_ebil ! PRINT level for energy conserv. diag. SAVE ip_ebil DATA ip_ebil/2/ INTEGER if_ebil ! level for energy conserv. dignostics SAVE if_ebil c+jld ec_conser REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique REAL ZRCPD c-jld ec_conser cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels REAL t2m(klon,nbsrf), q2m(klon,nbsrf) !temperature, humidite a 2m REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max c c Declaration des constantes et des fonctions thermodynamiques c #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" c====================================================================== modname = 'physiq' IF (if_ebil.ge.1) THEN DO i=1,klon zero_v(i)=0. END DO END IF ok_sync=.TRUE. IF (nqmax .LT. 2) THEN PRINT*, 'eaux vapeur et liquide sont indispensables' CALL ABORT ENDIF IF (debut) THEN CALL suphec ! initialiser constantes et parametres phys. ENDIF c====================================================================== xjour = rjourvrai c c Si c'est le debut, il faut initialiser plusieurs choses c ******** c IF (debut) THEN C IF (if_ebil.ge.1) d_h_vcol_phy=0. c c appel a la lecture du run.def physique c call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, . ok_instan, fact_cldcon, facttemps,ok_newmicro, . iflag_cldcon,ratqsbas,ratqshaut, if_ebil) cIM . , RI0) c c c Initialiser les compteurs: c frugs = 0. itap = 0 itaprad = 0 c CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0, . rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow, . falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown, . dlw,radsol,frugs,agesno,clesphy0, . zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0, . t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon ) c radpas = NINT( 86400./dtime/nbapp_rad) c C on remet le calendrier à zero c IF (raz_date .eq. 1) THEN itau_phy = 0 ENDIF c CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe, , ok_instan, ok_region ) c IF (ABS(dtime-pdtphys).GT.0.001) THEN PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys abort_message=' See above ' call abort_gcm(modname,abort_message,1) ENDIF IF (nlon .NE. klon) THEN PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon abort_message=' See above ' call abort_gcm(modname,abort_message,1) ENDIF IF (nlev .NE. klev) THEN PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev abort_message=' See above ' call abort_gcm(modname,abort_message,1) ENDIF c IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN PRINT*, 'Nbre d appels au rayonnement insuffisant' PRINT*, "Au minimum 4 appels par jour si cycle diurne" abort_message=' See above ' call abort_gcm(modname,abort_message,1) ENDIF PRINT*, "Clef pour la convection, iflag_con=", iflag_con PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl c cKE43 c Initialisation pour la convection de K.E. (sb): IF (iflag_con.GE.3) THEN PRINT*, "*** Convection de Kerry Emanuel 4.3 " PRINT*, "On va utiliser le melange convectif des traceurs qui" PRINT*, "est calcule dans convect4.3" PRINT*, " !!! penser aux logical flags de phytrac" DO i = 1, klon ema_cbmf(i) = 0. ema_pcb(i) = 0. ema_pct(i) = 0. ema_workcbmf(i) = 0. ENDDO cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG DO i = 1, klon ibas_con(i) = 1 itop_con(i) = klev+1 ENDDO cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END ENDIF c34EK IF (ok_orodr) THEN DO i=1,klon rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) ENDDO CALL SUGWD(klon,klev,paprs,pplay) DO i=1,klon zuthe(i)=0. zvthe(i)=0. if(zstd(i).gt.10.)then zuthe(i)=(1.-zgam(i))*cos(zthe(i)) zvthe(i)=(1.-zgam(i))*sin(zthe(i)) endif ENDDO ENDIF c c lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours PRINT*,'La frequence de lecture surface est de ', lmt_pas c ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours IF (ok_mensuel) THEN PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth ENDIF ecrit_day = NINT(86400./dtime *1.0) ! tous les jours IF (ok_journe) THEN PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day ENDIF ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures IF (ok_instan) THEN PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins ENDIF ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour IF (ok_region) THEN PRINT*, 'La frequence de sortie region est de ', ecrit_reg ENDIF c c Initialiser le couplage si necessaire c npas = 0 nexca = 0 if (ocean == 'couple') then npas = itaufin/ iphysiq nexca = 86400 / dtime write(*,*)' ##### Ocean couple #####' write(*,*)' Valeurs des pas de temps' write(*,*)' npas = ', npas write(*,*)' nexca = ', nexca endif c c cccIM capemaxcels = 't_max(X)' t2mincels = 't_min(X)' t2maxcels = 't_max(X)' cccIM cf. FH c c============================================================= c Initialisation des sorties c============================================================= #ifdef histhf #include "ini_histhf.h" #endif #include "ini_histday.h" #include "ini_histmth.h" #undef histmthNMC #define histmthNMC #ifdef histmthNMC #include "ini_histmthNMC.h" #endif #include "ini_histins.h" #ifdef histREGDYN #include "ini_histREGDYN.h" #endif #ifdef histISCCP #include "ini_histISCCP.h" #endif cXXXPB Positionner date0 pour initialisation de ORCHIDEE date0 = zjulian C date0 = day_ini WRITE(*,*) 'physiq date0 : ',date0 c c c c Prescrire l'ozone dans l'atmosphere c c cc DO i = 1, klon cc DO k = 1, klev cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20) cc ENDDO cc ENDDO c c ENDIF c c **************** Fin de IF ( debut ) *************** c c c Mettre a zero des variables de sortie (pour securite) c DO i = 1, klon d_ps(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon d_t(i,k) = 0.0 d_u(i,k) = 0.0 d_v(i,k) = 0.0 ENDDO ENDDO DO iq = 1, nqmax DO k = 1, klev DO i = 1, klon d_qx(i,k,iq) = 0.0 ENDDO ENDDO ENDDO c c Ne pas affecter les valeurs entrees de u, v, h, et q c DO k = 1, klev DO i = 1, klon t_seri(i,k) = t(i,k) u_seri(i,k) = u(i,k) v_seri(i,k) = v(i,k) q_seri(i,k) = qx(i,k,ivap) ql_seri(i,k) = qx(i,k,iliq) qs_seri(i,k) = 0. ENDDO ENDDO IF (nqmax.GE.3) THEN DO iq = 3, nqmax DO k = 1, klev DO i = 1, klon tr_seri(i,k,iq-2) = qx(i,k,iq) ENDDO ENDDO ENDDO ELSE DO k = 1, klev DO i = 1, klon tr_seri(i,k,1) = 0.0 ENDDO ENDDO ENDIF C DO i = 1, klon ztsol(i) = 0. ENDDO DO nsrf = 1, nbsrf DO i = 1, klon ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO C IF (if_ebil.ge.1) THEN ztit='after dynamic' CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) C Comme les tendances de la physique sont ajoute dans la dynamique, C on devrait avoir que la variation d'entalpie par la dynamique C est egale a la variation de la physique au pas de temps precedent. C Donc la somme de ces 2 variations devrait etre nulle. call diagphy(paire,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, zero_v e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol+d_h_vcol_phy, d_qt, 0. s , fs_bound, fq_bound ) END IF c Diagnostiquer la tendance dynamique c IF (ancien_ok) THEN DO k = 1, klev DO i = 1, klon d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime ENDDO ENDDO ELSE DO k = 1, klev DO i = 1, klon d_t_dyn(i,k) = 0.0 d_q_dyn(i,k) = 0.0 ENDDO ENDDO ancien_ok = .TRUE. ENDIF c c Ajouter le geopotentiel du sol: c DO k = 1, klev DO i = 1, klon zphi(i,k) = pphi(i,k) + pphis(i) ENDDO ENDDO c c Verifier les temperatures c CALL hgardfou(t_seri,ftsol,'debutphy') c c Incrementer le compteur de la physique c itap = itap + 1 julien = MOD(NINT(xjour),360) if (julien .eq. 0) julien = 360 c c Mettre en action les conditions aux limites (albedo, sst, etc.). c Prescrire l'ozone et calculer l'albedo sur l'ocean. c IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN PRINT *,' PHYS cond julien ',julien CALL ozonecm( FLOAT(julien), rlat, paprs, wo) ENDIF c c Re-evaporer l'eau liquide nuageuse c DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse DO i = 1, klon zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k)) c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k)) zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k)) zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k))) zb = MAX(0.0,ql_seri(i,k)) za = - MAX(0.0,ql_seri(i,k)) . * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta) t_seri(i,k) = t_seri(i,k) + za q_seri(i,k) = q_seri(i,k) + zb ql_seri(i,k) = 0.0 d_t_eva(i,k) = za d_q_eva(i,k) = zb ENDDO ENDDO c IF (if_ebil.ge.2) THEN ztit='after reevap' CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(paire,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, zero_v e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) C END IF C c c Appeler la diffusion verticale (programme de couche limite) c DO i = 1, klon c if (.not. ok_veget) then c frugs(i,is_ter) = SQRT(frugs(i,is_ter)**2+rugoro(i)**2) c endif c frugs(i,is_lic) = rugoro(i) c frugs(i,is_oce) = rugmer(i) c frugs(i,is_sic) = 0.001 zxrugs(i) = 0.0 ENDDO DO nsrf = 1, nbsrf DO i = 1, klon c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001) frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015) ENDDO ENDDO DO nsrf = 1, nbsrf DO i = 1, klon zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO c C calculs necessaires au calcul de l'albedo dans l'interface c CALL orbite(FLOAT(julien),zlongi,dist) IF (cycle_diurne) THEN zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s) CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract) ELSE rmu0 = -999.999 ENDIF cIM BEG DO i=1, klon sunlit(i)=1 IF(rmu0(i).EQ.0.) sunlit(i)=0 nbsunlit(1,i)=FLOAT(sunlit(i)) ENDDO cIM END C Calcul de l'abedo moyen par maille albsol(:)=0. albsollw(:)=0. DO nsrf = 1, nbsrf DO i = 1, klon albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf) albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf) ENDDO ENDDO C C Repartition sous maille des flux LW et SW C Modif OM+PASB+JLD C Repartition du longwave par sous-surface linearisee Cn DO nsrf = 1, nbsrf DO i = 1, klon c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4 c$$$ fsollw(i,nsrf) = sollw(i) fsollw(i,nsrf) = sollw(i) $ + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf)) fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i)) ENDDO ENDDO fder = dlw CALL clmain(dtime,itap,date0,pctsrf, e t_seri,q_seri,u_seri,v_seri, e julien, rmu0, co2_ppm, e ok_veget, ocean, npas, nexca, ftsol, $ soil_model,cdmmax, cdhmax, $ ksta, ksta_ter, ok_kzmin, ftsoil, qsol, $ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw, $ fluxlat, cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder, e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder, e rlon, rlat, cufi, cvfi, frugs, e debut, lafin, agesno,rugoro , s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts, s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm, s dsens, devap, cIM cf JLD s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m) s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m, s fqcalving,ffonte) c CXXX PB CXXX Incrementation des flux CXXX zxfluxt=0. zxfluxq=0. zxfluxu=0. zxfluxv=0. DO nsrf = 1, nbsrf DO k = 1, klev DO i = 1, klon zxfluxt(i,k) = zxfluxt(i,k) + $ fluxt(i,k,nsrf) * pctsrf( i, nsrf) zxfluxq(i,k) = zxfluxq(i,k) + $ fluxq(i,k,nsrf) * pctsrf( i, nsrf) zxfluxu(i,k) = zxfluxu(i,k) + $ fluxu(i,k,nsrf) * pctsrf( i, nsrf) zxfluxv(i,k) = zxfluxv(i,k) + $ fluxv(i,k,nsrf) * pctsrf( i, nsrf) END DO END DO END DO DO i = 1, klon sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol fder(i) = dlw(i) + dsens(i) + devap(i) ENDDO DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k) q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k) u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k) v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k) ENDDO ENDDO c IF (if_ebil.ge.2) THEN ztit='after clmain' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(paire,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, sens e , evap , zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF C c c Incrementer la temperature du sol c DO i = 1, klon zxtsol(i) = 0.0 zxfluxlat(i) = 0.0 cccIM zt2m(i) = 0.0 zq2m(i) = 0.0 zu10m(i) = 0.0 zv10m(i) = 0.0 cIM cf JLD ?? zxffonte(i) = 0.0 zxfqcalving(i) = 0.0 c IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + $ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA) $ THEN WRITE(*,*) 'physiq : pb sous surface au point ', i, $ pctsrf(i, 1 : nbsrf) ENDIF ENDDO DO nsrf = 1, nbsrf DO i = 1, klon c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf) cIM cf. JLD wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf) $ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf) zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf) zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf) cccIM zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf) zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf) zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf) zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf) cIM cf JLD ?? zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf) zxfqcalving(i) = zxfqcalving(i) + . fqcalving(i,nsrf)*pctsrf(i,nsrf) c ENDIF ENDDO ENDDO c c Si une sous-fraction n'existe pas, elle prend la temp. moyenne c DO nsrf = 1, nbsrf DO i = 1, klon IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i) cccIM IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i) IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i) IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i) IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i) cIM cf JLD ?? IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i) IF (pctsrf(i,nsrf) .LT. epsfra) . fqcalving(i,nsrf) = zxfqcalving(i) ENDDO ENDDO c c c Calculer la derive du flux infrarouge c cXXX DO nsrf = 1, nbsrf DO i = 1, klon cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3 cXXX . *(ftsol(i,nsrf)-zxtsol(i)) cXXX . *pctsrf(i,nsrf) cXXX ENDIF cXXX ENDDO ENDDO c c Appeler la convection (au choix) c DO k = 1, klev DO i = 1, klon conv_q(i,k) = d_q_dyn(i,k) . + d_q_vdf(i,k)/dtime conv_t(i,k) = d_t_dyn(i,k) . + d_t_vdf(i,k)/dtime ENDDO ENDDO IF (check) THEN za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire) PRINT*, "avantcon=", za ENDIF zx_ajustq = .FALSE. IF (iflag_con.EQ.2) zx_ajustq=.TRUE. IF (zx_ajustq) THEN DO i = 1, klon z_avant(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) . *(paprs(i,k)-paprs(i,k+1))/RG ENDDO ENDDO ENDIF IF (iflag_con.EQ.1) THEN stop'reactiver le call conlmd dans physiq.F' c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q, c . d_t_con, d_q_con, c . rain_con, snow_con, ibas_con, itop_con) ELSE IF (iflag_con.EQ.2) THEN CALL conflx(dtime, paprs, pplay, t_seri, q_seri, e conv_t, conv_q, zxfluxq(1,1), omega, s d_t_con, d_q_con, rain_con, snow_con, s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, s kcbot, kctop, kdtop, pmflxr, pmflxs) WHERE (rain_con < 0.) rain_con = 0. WHERE (snow_con < 0.) snow_con = 0. DO i = 1, klon ibas_con(i) = klev+1 - kcbot(i) itop_con(i) = klev+1 - kctop(i) ENDDO ELSE IF (iflag_con.GE.3) THEN c nb of tracers for the KE convection: if (nqmax .GE. 4) then ntra = nbtr else ntra = 1 endif c c sb, oct02: c Schema de convection modularise et vectorise: c (driver commun aux versions 3 et 4) c IF (ok_cvl) THEN ! new driver for convectL CALL concvl (iflag_con, . dtime,paprs,pplay,t_seri,q_seri, . u_seri,v_seri,tr_seri,nbtr, . ema_work1,ema_work2, . d_t_con,d_q_con,d_u_con,d_v_con,d_tr, . rain_con, snow_con, ibas_con, itop_con, . upwd,dnwd,dnwd0, . Ma,cape,tvp,iflagctrl, . pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd) cIM cf. FH clwcon0=qcondc ELSE ! ok_cvl c print*,'Avant conema OUI' CALL conema3 (dtime, . paprs,pplay,t_seri,q_seri, . u_seri,v_seri,tr_seri,nbtr, . ema_work1,ema_work2, . d_t_con,d_q_con,d_u_con,d_v_con,d_tr, . rain_con, snow_con, ibas_con, itop_con, . upwd,dnwd,dnwd0,bas,top, . Ma,cape,tvp,rflag, . pbase . ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr . ,clwcon0) print*,'Apres conema3 ' ENDIF ! ok_cvl IF (.NOT. ok_gust) THEN do i = 1, klon wd(i)=0.0 enddo ENDIF c =================================================================== c c Calcul des proprietes des nuages convectifs c DO k = 1, klev DO i = 1, klon zx_t = t_seri(i,k) IF (thermcep) THEN zdelta = MAX(0.,SIGN(1.,rtt-zx_t)) zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k) zx_qs = MIN(0.5,zx_qs) zcor = 1./(1.-retv*zx_qs) zx_qs = zx_qs*zcor ELSE IF (zx_t.LT.t_coup) THEN zx_qs = qsats(zx_t)/pplay(i,k) ELSE zx_qs = qsatl(zx_t)/pplay(i,k) ENDIF ENDIF zqsat(i,k)=zx_qs ENDDO ENDDO c calcul des proprietes des nuages convectifs clwcon0(:,:)=fact_cldcon*clwcon0(:,:) call clouds_gno s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0) c =================================================================== c DO i = 1, klon ema_pcb(i) = pbase(i) ENDDO DO i = 1, klon ema_pct(i) = paprs(i,itop_con(i)) ENDDO DO i = 1, klon ema_cbmf(i) = ema_workcbmf(i) ENDDO ELSE PRINT*, "iflag_con non-prevu", iflag_con CALL abort ENDIF c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri, c . d_u_con, d_v_con) DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_con(i,k) q_seri(i,k) = q_seri(i,k) + d_q_con(i,k) u_seri(i,k) = u_seri(i,k) + d_u_con(i,k) v_seri(i,k) = v_seri(i,k) + d_v_con(i,k) ENDDO ENDDO c IF (if_ebil.ge.2) THEN ztit='after convect' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(paire,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, zero_v e , zero_v, rain_con, snow_con, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF C IF (check) THEN za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire) PRINT*, "aprescon=", za zx_t = 0.0 za = 0.0 DO i = 1, klon za = za + paire(i)/FLOAT(klon) zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon) ENDDO zx_t = zx_t/za*dtime PRINT*, "Precip=", zx_t ENDIF IF (zx_ajustq) THEN DO i = 1, klon z_apres(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) . *(paprs(i,k)-paprs(i,k+1))/RG ENDDO ENDDO DO i = 1, klon z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) . /z_apres(i) ENDDO DO k = 1, klev DO i = 1, klon IF (z_factor(i).GT.(1.0+1.0E-08) .OR. . z_factor(i).LT.(1.0-1.0E-08)) THEN q_seri(i,k) = q_seri(i,k) * z_factor(i) ENDIF ENDDO ENDDO ENDIF zx_ajustq=.FALSE. c IF (nqmax.GT.2) THEN !--melange convectif de traceurs c IF (iflag_con .NE. 2 .AND. debut) THEN PRINT*, 'Pour l instant, seul conflx fonctionne ', $ 'avec traceurs', iflag_con PRINT*,' Mettre iflag_con', $ ' = 2 dans run.def et repasser' c CALL abort ENDIF c ENDIF !--nqmax.GT.2 c c Appeler l'ajustement sec c CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs) DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k) q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k) ENDDO ENDDO c IF (if_ebil.ge.2) THEN ztit='after dry_adjust' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) END IF c------------------------------------------------------------------------- c Caclul des ratqs c------------------------------------------------------------------------- c print*,'calcul des ratqs' c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q c ---------------- c on ecrase le tableau ratqsc calcule par clouds_gno if (iflag_cldcon.eq.1) then do k=1,klev do i=1,klon if(ptconv(i,k)) then ratqsc(i,k)=ratqsbas s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k) else ratqsc(i,k)=0. endif enddo enddo endif c ratqs stables c ------------- do k=1,klev ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)* s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.) enddo c ratqs final c ----------- if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then c les ratqs sont une conbinaison de ratqss et ratqsc c ratqs final c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de c relaxation des ratqs c facttemps=exp(-pdtphys/1.e4) facteur=exp(-pdtphys*facttemps) ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:)) ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:)) c print*,'calcul des ratqs fini' else c on ne prend que le ratqs stable pour fisrtilp ratqs(:,:)=ratqss(:,:) endif c c Appeler le processus de condensation a grande echelle c et le processus de precipitation c------------------------------------------------------------------------- CALL fisrtilp(dtime,paprs,pplay, . t_seri, q_seri,ptconv,ratqs, . d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, . rain_lsc, snow_lsc, . pfrac_impa, pfrac_nucl, pfrac_1nucl, . frac_impa, frac_nucl, . prfl, psfl, rhcl) WHERE (rain_lsc < 0) rain_lsc = 0. WHERE (snow_lsc < 0) snow_lsc = 0. DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k) q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k) ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k) cldfra(i,k) = rneb(i,k) IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k) ENDDO ENDDO IF (check) THEN za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire) PRINT*, "apresilp=", za zx_t = 0.0 za = 0.0 DO i = 1, klon za = za + paire(i)/FLOAT(klon) zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon) ENDDO zx_t = zx_t/za*dtime PRINT*, "Precip=", zx_t ENDIF c IF (if_ebil.ge.2) THEN ztit='after fisrt' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(paire,ztit,ip_ebil e , zero_v, zero_v, zero_v, zero_v, zero_v e , zero_v, rain_lsc, snow_lsc, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF c c------------------------------------------------------------------- c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT c------------------------------------------------------------------- c 1. NUAGES CONVECTIFS c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke c Nuages diagnostiques pour Tiedtke CALL diagcld1(paprs,pplay, . rain_con,snow_con,ibas_con,itop_con, . diafra,dialiq) DO k = 1, klev DO i = 1, klon IF (diafra(i,k).GT.cldfra(i,k)) THEN cldliq(i,k) = dialiq(i,k) cldfra(i,k) = diafra(i,k) ENDIF ENDDO ENDDO ELSE IF (iflag_cldcon.eq.3) THEN c On prend pour les nuages convectifs le max du calcul de la c convection et du calcul du pas de temps précédent diminué d'un facteur c facttemps c facttemps=pdtphys/1.e4 facteur = pdtphys *facttemps do k=1,klev do i=1,klon rnebcon(i,k)=rnebcon(i,k)*facteur if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k)) s then rnebcon(i,k)=rnebcon0(i,k) clwcon(i,k)=clwcon0(i,k) endif enddo enddo cIM calcul nuages par le simulateur ISCCP IF (ok_isccp) THEN cIM calcul tau. emi nuages convectifs convfra(:,:)=rnebcon(:,:) convliq(:,:)=rnebcon(:,:)*clwcon(:,:) CALL newmicro (paprs, pplay,ok_newmicro, . t_seri, convliq, convfra, dtau_c, dem_c, . cldh_c, cldl_c, cldm_c, cldt_c, cldq_c, . flwp_c, fiwp_c, flwc_c, fiwc_c) c cIM calcul tau. emi nuages startiformes CALL newmicro (paprs, pplay,ok_newmicro, . t_seri, cldliq, cldfra, dtau_s, dem_s, . cldh_s, cldl_s, cldm_s, cldt_s, cldq_s, . flwp_s, fiwp_s, flwc_s, fiwc_s) c cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.) cIM inversion des niveaux de pression ==> de haut en bas CALL haut2bas(klon, klev, pplay, pfull) CALL haut2bas(klon, klev, q_seri, qv) CALL haut2bas(klon, klev, cldtot, cc) CALL haut2bas(klon, klev, rnebcon, conv) CALL haut2bas(klon, klev, dtau_s, dtau_sH2B) CALL haut2bas(klon, klev, dtau_c, dtau_cH2B) CALL haut2bas(klon, klev, t_seri, at) CALL haut2bas(klon, klev, dem_s, dem_sH2B) CALL haut2bas(klon, klev, dem_c, dem_cH2B) CALL haut2bas(klon, klevp1, paprs, phalf) c open(99,file='tautab.bin',access='sequential', c $ form='unformatted',status='old') c read(99) tautab cIM210503 IF (debut) THEN open(99,file='tautab.formatted', FORM='FORMATTED') read(99,'(f30.20)') tautab close(99) c open(99,file='invtau.formatted',form='FORMATTED') read(99,'(i10)') invtau close(99) c cIM: calcul coordonnees regions pour statistiques distribution cIM: nuages en ftion du regime dynamique pour regions oceaniques IF (ok_regdyn) THEN !histREGDYN nsrf=3 DO nreg=1, nbregdyn DO i=1, klon c IF (debut) THEN IF(rlon(i).LT.0.) THEN rlonPOS(i)=rlon(i)+360. ELSE rlonPOS(i)=rlon(i) ENDIF c ENDIF pct_ocean(i,nreg)=0 c test si c'est 1 point d'ocean IF(pctsrf(i,nsrf).EQ.1.) THEN IF(nreg.EQ.1) THEN c TROP IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN pct_ocean(i,nreg)=1 ENDIF c PACIFIQUE NORD ELSEIF(nreg.EQ.2) THEN IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN pct_ocean(i,nreg)=1 ENDIF ENDIF c CALIFORNIE ST-CU ELSEIF(nreg.EQ.3) THEN IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN pct_ocean(i,nreg)=1 ENDIF ENDIF c HAWAI ELSEIF(nreg.EQ.4) THEN IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN pct_ocean(i,nreg)=1 ENDIF ENDIF c WARM POOL ELSEIF(nreg.EQ.5) THEN IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN pct_ocean(i,nreg)=1 ENDIF ENDIF ENDIF !nbregdyn c TROP c IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN c pct_ocean(i)=.TRUE. c WRITE(*,*) 'pct_ocean =',i, rlon(i), rlat(i) c ENDIF !lon c ENDIF !lat ENDIF !pctsrf ENDDO !klon ENDDO !nbregdyn ENDIF !ok_regdyn cIM somme de toutes les nhistoW BEG DO nreg = 1, nbregdyn DO k = 1, kmaxm1 DO l = 1, lmaxm1 DO iw = 1, iwmax nhistoWt(k,l,iw,nreg)=0. ENDDO !iw ENDDO !l ENDDO !k ENDDO !nreg cIM somme de toutes les nhistoW END c cIM: initialisation de seed DO i=1, klon seed(i)=i+100 ENDDO ENDIF !debut cIM: pas de debug, debugcol debug=0 debugcol=0 cIM260503 c o500 ==> distribution nuage ftion du regime dynamique a 500 hPa DO k=1, klevm1 kp1=k+1 c PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1) if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN DO i=1, klon o500(i)=omega(i,k)*RDAY/100. c if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1) ENDDO GOTO 1000 endif 1000 continue ENDDO CALL ISCCP_CLOUD_TYPES( & debug, & debugcol, & klon, & sunlit, & klev, & ncol, & seed, & pfull, & phalf, & qv, cc, conv, dtau_sH2B, dtau_cH2B, & top_height, & overlap, & tautab, & invtau, & ztsol, & emsfc_lw, & at, dem_sH2B, dem_cH2B, & fq_isccp, & totalcldarea, & meanptop, & meantaucld, & boxtau, & boxptop) c passage de la grille (klon,7,7) a (iim,jjmp1,7,7) DO l=1, lmaxm1 DO k=1, kmaxm1 DO i=1, iim fq4d(i,1,k,l)=fq_isccp(1,k,l) ENDDO DO j=2, jjm DO i=1, iim ig=i+1+(j-2)*iim fq4d(i,j,k,l)=fq_isccp(ig,k,l) ENDDO ENDDO DO i=1, iim fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l) ENDDO ENDDO ENDDO c DO l=1, lmaxm1 DO k=1, kmaxm1 DO j=1, jjmp1 DO i=1, iim ni=(i-1)*lmaxm1+l nj=(j-1)*kmaxm1+k fq3d(ni,nj)=fq4d(i,j,k,l) ENDDO ENDDO ENDDO ENDDO c c calculs statistiques distribution nuage ftion du regime dynamique c c Ce calcul doit etre fait a partir de valeurs mensuelles ?? CALL histo_o500_pctau(nbregdyn,pct_ocean,o500,fq_isccp, &histoW,nhistoW) c c nhistoWt = somme de toutes les nhistoW DO nreg=1, nbregdyn DO k = 1, kmaxm1 DO l = 1, lmaxm1 DO iw = 1, iwmax nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+ & nhistoW(k,l,iw,nreg) ENDDO ENDDO ENDDO ENDDO c ENDIF !ok_isccp c On prend la somme des fractions nuageuses et des contenus en eau cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.) cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:) ENDIF c c 2. NUAGES STARTIFORMES c IF (ok_stratus) THEN CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq) DO k = 1, klev DO i = 1, klon IF (diafra(i,k).GT.cldfra(i,k)) THEN cldliq(i,k) = dialiq(i,k) cldfra(i,k) = diafra(i,k) ENDIF ENDDO ENDDO ENDIF c c Precipitation totale c DO i = 1, klon rain_fall(i) = rain_con(i) + rain_lsc(i) snow_fall(i) = snow_con(i) + snow_lsc(i) ENDDO c IF (if_ebil.ge.2) THEN ztit="after diagcld" CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) END IF c c Calculer l'humidite relative pour diagnostique c DO k = 1, klev DO i = 1, klon zx_t = t_seri(i,k) IF (thermcep) THEN zdelta = MAX(0.,SIGN(1.,rtt-zx_t)) zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k) zx_qs = MIN(0.5,zx_qs) zcor = 1./(1.-retv*zx_qs) zx_qs = zx_qs*zcor ELSE IF (zx_t.LT.t_coup) THEN zx_qs = qsats(zx_t)/pplay(i,k) ELSE zx_qs = qsatl(zx_t)/pplay(i,k) ENDIF ENDIF zx_rh(i,k) = q_seri(i,k)/zx_qs zqsat(i,k)=zx_qs ENDDO ENDDO c c Calculer les parametres optiques des nuages et quelques c parametres pour diagnostiques: c if (ok_newmicro) then CALL newmicro (paprs, pplay,ok_newmicro, . t_seri, cldliq, cldfra, cldtau, cldemi, . cldh, cldl, cldm, cldt, cldq, . flwp, fiwp, flwc, fiwc) else CALL nuage (paprs, pplay, . t_seri, cldliq, cldfra, cldtau, cldemi, . cldh, cldl, cldm, cldt, cldq) endif c c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol. c IF (MOD(itaprad,radpas).EQ.0) THEN DO i = 1, klon albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce) . + falbe(i,is_lic) * pctsrf(i,is_lic) . + falbe(i,is_ter) * pctsrf(i,is_ter) . + falbe(i,is_sic) * pctsrf(i,is_sic) albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce) . + falblw(i,is_lic) * pctsrf(i,is_lic) . + falblw(i,is_ter) * pctsrf(i,is_ter) . + falblw(i,is_sic) * pctsrf(i,is_sic) ENDDO ! if (debut) then ! albsol1 = albsol ! albsollw1 = albsollw ! endif ! albsol = albsol1 ! albsollw = albsollw1 CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS) e (dist, rmu0, fract, e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri, e wo, e cldfra, cldemi, cldtau, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, s sollwdown, s topsw0,toplw0,solsw0,sollw0, s swdn0, swdn, swup0, swup ) itaprad = 0 ENDIF itaprad = itaprad + 1 c c Ajouter la tendance des rayonnements (tous les pas) c DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) . + (heat(i,k)-cool(i,k)) * dtime/86400. ENDDO ENDDO c IF (if_ebil.ge.2) THEN ztit='after rad' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) call diagphy(paire,ztit,ip_ebil e , topsw, toplw, solsw, sollw, zero_v e , zero_v, zero_v, zero_v, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) END IF c c c Calculer l'hydrologie de la surface c c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap, c . agesno, ftsol,fqsurf,fsnow, ruis) c DO i = 1, klon zxqsurf(i) = 0.0 zxsnow(i) = 0.0 ENDDO DO nsrf = 1, nbsrf DO i = 1, klon zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf) zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf) ENDDO ENDDO c c Si une sous-fraction n'existe pas, elle prend la valeur moyenne c cXXX DO nsrf = 1, nbsrf cXXX DO i = 1, klon cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN cXXX fqsurf(i,nsrf) = zxqsurf(i) cXXX fsnow(i,nsrf) = zxsnow(i) cXXX ENDIF cXXX ENDDO cXXX ENDDO c c Calculer le bilan du sol et la derive de temperature (couplage) c DO i = 1, klon c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT c a la demande de JLD bils(i) = radsol(i) - sens(i) + zxfluxlat(i) ENDDO c cmoddeblott(jan95) c Appeler le programme de parametrisation de l'orographie c a l'echelle sous-maille: c IF (ok_orodr) THEN c c selection des points pour lesquels le shema est actif: igwd=0 DO i=1,klon itest(i)=0 c IF ((zstd(i).gt.10.0)) THEN IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN itest(i)=1 igwd=igwd+1 idx(igwd)=i ENDIF ENDDO c igwdim=MAX(1,igwd) c CALL drag_noro(klon,klev,dtime,paprs,pplay, e zmea,zstd, zsig, zgam, zthe,zpic,zval, e igwd,idx,itest, e t_seri, u_seri, v_seri, s zulow, zvlow, zustr, zvstr, s d_t_oro, d_u_oro, d_v_oro) c c ajout des tendances DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k) u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k) v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k) ENDDO ENDDO c ENDIF ! fin de test sur ok_orodr c IF (ok_orolf) THEN c c selection des points pour lesquels le shema est actif: igwd=0 DO i=1,klon itest(i)=0 IF ((zpic(i)-zmea(i)).GT.100.) THEN itest(i)=1 igwd=igwd+1 idx(igwd)=i ENDIF ENDDO c igwdim=MAX(1,igwd) c CALL lift_noro(klon,klev,dtime,paprs,pplay, e rlat,zmea,zstd,zpic, e itest, e t_seri, u_seri, v_seri, s zulow, zvlow, zustr, zvstr, s d_t_lif, d_u_lif, d_v_lif) c c ajout des tendances DO k = 1, klev DO i = 1, klon t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k) u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k) v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k) ENDDO ENDDO c ENDIF ! fin de test sur ok_orolf c IF (if_ebil.ge.2) THEN ztit='after orography' CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) END IF c c cAA cAA Installation de l'interface online-offline pour traceurs cAA c==================================================================== c Calcul des tendances traceurs c==================================================================== C Pascale : il faut quand meme apeller phytrac car il gere les sorties cKE43 des traceurs => il faut donc mettre des flags a .false. IF (iflag_con.GE.3) THEN c on ajoute les tendances calculees par KE43 cXXX OM on onhibe la convection sur les traceurs DO iq=1, nqmax-2 ! Sandrine a -3 ??? cXXX OM on inhibe la convection sur les traceur cXXX DO k = 1, nlev cXXX DO i = 1, klon cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq) cXXX ENDDO cXXX ENDDO WRITE(iqn,'(i2.2)') iq CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn) ENDDO CMAF modif pour garder info du nombre de traceurs auxquels C la physique s'applique ELSE CMAF modif pour garder info du nombre de traceurs auxquels C la physique s'applique C call phytrac (rnpb, I debut,lafin, I nqmax-2, I nlon,nlev,dtime, I t,paprs,pplay, I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, I ycoefh,yu1,yv1,ftsol,pctsrf,rlat, I frac_impa, frac_nucl, I rlon,presnivs,paire,pphis, O tr_seri) ENDIF IF (offline) THEN call phystokenc ( I nlon,nlev,pdtphys,rlon,rlat, I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, I ycoefh,yu1,yv1,ftsol,pctsrf, I frac_impa, frac_nucl, I pphis,paire,dtime,itap) ENDIF c c Calculer le transport de l'eau et de l'energie (diagnostique) c CALL transp (paprs,zxtsol, e t_seri, q_seri, u_seri, v_seri, zphi, s ve, vq, ue, uq) c c c Accumuler les variables a stocker dans les fichiers histoire: c c c c+jld ec_conser DO k = 1, klev DO i = 1, klon ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i,k)) d_t_ec(i,k)=0.5/ZRCPD $ *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2) t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k) d_t_ec(i,k) = d_t_ec(i,k)/dtime END DO END DO c-jld ec_conser IF (if_ebil.ge.1) THEN ztit='after physic' CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec) C Comme les tendances de la physique sont ajoute dans la dynamique, C on devrait avoir que la variation d'entalpie par la dynamique C est egale a la variation de la physique au pas de temps precedent. C Donc la somme de ces 2 variations devrait etre nulle. call diagphy(paire,ztit,ip_ebil e , topsw, toplw, solsw, sollw, sens e , evap, rain_fall, snow_fall, ztsol e , d_h_vcol, d_qt, d_ec s , fs_bound, fq_bound ) C d_h_vcol_phy=d_h_vcol C END IF C c======================================================================= c SORTIES c======================================================================= c Interpollation sur quelques niveaux de pression c ----------------------------------------------- c cIM sorties sur les 17 niveaux de pression du NMC c 1000 hPa call plevel(klon,klev,.true. ,pplay,100000.,u_seri,u1000) call plevel(klon,klev,.false.,pplay,100000.,v_seri,v1000) c 925 hPa call plevel(klon,klev,.true. ,pplay,92500.,u_seri,u925) call plevel(klon,klev,.false.,pplay,92500.,v_seri,v925) c 850 hPa call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850) call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850) c 700 hPa call plevel(klon,klev,.true. ,pplay,70000.,u_seri,u700) call plevel(klon,klev,.false.,pplay,70000.,v_seri,v700) c 600 hPa call plevel(klon,klev,.true. ,pplay,60000.,u_seri,u600) call plevel(klon,klev,.false.,pplay,60000.,v_seri,v600) c 500 hPa call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500) call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500) c 400 hPa call plevel(klon,klev,.true. ,pplay,40000.,u_seri,u400) call plevel(klon,klev,.false.,pplay,40000.,v_seri,v400) c 300 hPa call plevel(klon,klev,.true. ,pplay,30000.,u_seri,u300) call plevel(klon,klev,.false.,pplay,30000.,v_seri,v300) c 250 hPa call plevel(klon,klev,.true. ,pplay,25000.,u_seri,u250) call plevel(klon,klev,.false.,pplay,25000.,v_seri,v250) c 200 hPa call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200) call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200) c 150 hPa call plevel(klon,klev,.true. ,pplay,15000.,u_seri,u150) call plevel(klon,klev,.false.,pplay,15000.,v_seri,v150) c 100 hPa call plevel(klon,klev,.true. ,pplay,10000.,u_seri,u100) call plevel(klon,klev,.false.,pplay,10000.,v_seri,v100) c 70 hPa call plevel(klon,klev,.true. ,pplay,7000.,u_seri,u70) call plevel(klon,klev,.false.,pplay,7000.,v_seri,v70) c 50 hPa call plevel(klon,klev,.true. ,pplay,5000.,u_seri,u50) call plevel(klon,klev,.false.,pplay,5000.,v_seri,v50) c 30 hPa call plevel(klon,klev,.true. ,pplay,3000.,u_seri,u30) call plevel(klon,klev,.false.,pplay,3000.,v_seri,v30) c 20 hPa call plevel(klon,klev,.true. ,pplay,2000.,u_seri,u20) call plevel(klon,klev,.false.,pplay,2000.,v_seri,v20) c 10 hPa call plevel(klon,klev,.true. ,pplay,1000.,u_seri,u10) call plevel(klon,klev,.false.,pplay,1000.,v_seri,v10) c call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500) call plevel(klon,klev,.true. ,paprs,50000.,omega,w500) c slp sea level pressure slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1))) c ccc prw = eau precipitable DO i = 1, klon prw(i) = 0. DO k = 1, klev prw(i) = prw(i) + . q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG ENDDO ENDDO c cIM sorties bilans energie cinetique et potentielle MJO DO k = 1, klev DO i = 1, klon d_u_oli(i,k) = d_u_oro(i,k) + d_u_lif(i,k) d_v_oli(i,k) = d_v_oro(i,k) + d_v_lif(i,k) ENDDO ENDDO c============================================================= c Ecriture des sorties c============================================================= #ifdef histREGDYN #include "write_histREGDYN.h" #endif #ifdef histISCCP #include "write_histISCCP.h" #endif #ifdef histhf #include "write_histhf.h" #endif #include "write_histday.h" #include "write_histmth.h" #ifdef histmthNMC #include "write_histmthNMC.h" #endif #include "write_histins.h" c============================================================= c c Convertir les incrementations en tendances c DO k = 1, klev DO i = 1, klon d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime ENDDO ENDDO c IF (nqmax.GE.3) THEN DO iq = 3, nqmax DO k = 1, klev DO i = 1, klon d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime ENDDO ENDDO ENDDO ENDIF c c Sauvegarder les valeurs de t et q a la fin de la physique: c DO k = 1, klev DO i = 1, klon t_ancien(i,k) = t_seri(i,k) q_ancien(i,k) = q_seri(i,k) ENDDO ENDDO c c==================================================================== c Si c'est la fin, il faut conserver l'etat de redemarrage c==================================================================== c IF (lafin) THEN itau_phy = itau_phy + itap ccc IF (ok_oasis) CALL quitcpl CALL phyredem ("restartphy.nc",dtime,radpas, . rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol, . fsnow, falbe,falblw, fevap, rain_fall, snow_fall, . solsw, sollwdown,dlw, . radsol,frugs,agesno, . zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro, . t_ancien, q_ancien, rnebcon, ratqs, clwcon) ENDIF RETURN END FUNCTION qcheck(klon,klev,paprs,q,ql,aire) IMPLICIT none c c Calculer et imprimer l'eau totale. A utiliser pour verifier c la conservation de l'eau c #include "YOMCST.h" INTEGER klon,klev REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev) REAL aire(klon) REAL qtotal, zx, qcheck INTEGER i, k c zx = 0.0 DO i = 1, klon zx = zx + aire(i) ENDDO qtotal = 0.0 DO k = 1, klev DO i = 1, klon qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) . *(paprs(i,k)-paprs(i,k+1))/RG ENDDO ENDDO c qcheck = qtotal/zx c RETURN END SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit) IMPLICIT none c c Tranformer une variable de la grille physique a c la grille d'ecriture c INTEGER nfield,nlon,iim,jjmp1, jjm REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield) c INTEGER i, n, ig c jjm = jjmp1 - 1 DO n = 1, nfield DO i=1,iim ecrit(i,n) = fi(1,n) ecrit(i+jjm*iim,n) = fi(nlon,n) ENDDO DO ig = 1, nlon - 2 ecrit(iim+ig,n) = fi(1+ig,n) ENDDO ENDDO RETURN END