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_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) cccIM INTEGER klevp1 PARAMETER(klevp1=klev+1) #include "raddim.h" cc REAL*8 ZFSUP(KDLON,KFLEV+1) cc REAL*8 ZFSDN(KDLON,KFLEV+1) cc REAL*8 ZFSUP0(KDLON,KFLEV+1) cc REAL*8 ZFSDN0(KDLON,KFLEV+1) c REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2) SAVE swdn0 , swdn, swup0, swup cccIM cf. FH real u850(klon),v850(klon),u200(klon),v200(klon) real u500(klon),v500(klon),phi500(klon),w500(klon) cIM real prw(klon) cIM ISCCP - proprietes microphysiques des nuages convectifs 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 kinv, linv cIM ISCCP simulator BEGIN INTEGER igfi2D(iim,jjmp1) cv3.4 INTEGER debug, debugcol INTEGER npoints PARAMETER(npoints=klon) INTEGER sunlit(klon) INTEGER ncol, seed(klon) cIM dans clesphys.h top_height, overlap c PARAMETER(ncol=100) c PARAMETER(ncol=625) PARAMETER(ncol=10) REAL tautab(0:255) INTEGER invtau(-20:45000) REAL emsfc_lw PARAMETER(emsfc_lw=0.99) REAL ran0 ! type for random number fuction REAL pfull(klon,klev) REAL phalf(klon,klev+1) REAL cldtot(klon,klev) REAL dtau_s(klon,klev) REAL dtau_c(klon,klev) REAL dem_s(klon,klev) REAL dem_c(klon,klev) cPLUS : variables de haut en bas pour le simulateur ISCCP 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 REAL fq_isccp(klon,7,7) REAL totalcldarea(klon) REAL meanptop(klon) REAL meantaucld(klon) REAL boxtau(klon,ncol) REAL boxptop(klon,ncol) c grille 4d physique INTEGER l, ni, nj, kmax, lmax, nrec INTEGER ni1, ni2, nj1, nj2 c PARAMETER(kmax=7, lmax=7) PARAMETER(kmax=8, lmax=8) INTEGER kmaxm1, lmaxm1 PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1) c INTEGER iimx7, jjmx7, jjmp1x7 c PARAMETER(iimx7=iim*7, jjmx7=jjm*7, jjmp1x7=jjmp1*7) c REAL fq4d(iim,jjmp1,7,7) c REAL fq3d(iimx7, jjmp1x7) INTEGER iimx8, jjmx8, jjmp1x8 PARAMETER(iimx8=iim*8, jjmx8=jjm*8, jjmp1x8=jjmp1*8) REAL fq4d(iim,jjmp1,8,8) REAL fq3d(iimx8, jjmp1x8) cIM180603 SAVE fq3d c REAL maxfq3d, minfq3d 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) INTEGER nreg, nbreg PARAMETER(nbreg=5) c REAL histoW(iwmax,kmaxm1,lmaxm1) REAL histoW(kmaxm1,lmaxm1,iwmax,nbreg) REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbreg) cIM180603 c SAVE histoW, nhistoW c SAVE nhistoW REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbreg) SAVE nhistoWt c REAL histoWinv(kmaxm1,lmaxm1,iwmax) c REAL nhistoW(kmaxm1,lmaxm1,iwmax) INTEGER linv c LOGICAL pct_ocean(klon,nbreg) INTEGER pct_ocean(klon,nbreg) REAL rlonPOS(klon) c CHARACTER*4 pdirect 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 REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax) c DATA zx_tau/0.1, 1.3, 3.6, 9.4, 23., 60./ c DATA zx_pc/50., 180., 310., 440., 560., 680., 800., 1015./ c DATA zx_pc/50., 180., 310., 440., 560., 680., 800./ cOK DATA zx_tau/0.0, 0.1, 1.3, 3.6, 9.4, 23., 60./ cOK DATA zx_pc/800., 680., 560., 440., 310., 180., 50./ c tester l'alure DATA zx_tau/1., 2., 3., 4., 5., 6., 7./ c DATA zx_pc/1., 2., 3., 4., 5., 6., 7./ DATA zx_pc/7., 6., 5., 4., 3., 2., 1./ INTEGER komega, nhoriRD c statistiques regime dynamique END c REAL del_lon(iim), del_lat(jjmp1) REAL del_lon, del_lat c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7) REAL zx_lonx8(iimx8), zx_latx8(jjmp1x8) c INTEGER nhorix7 INTEGER nhorix8 cIM ISCCP simulator END 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 cIM cf JLD 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 cIM cf. JLD 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 albsno ! calculer albedo sur neige 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 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 cccIM CHARACTER*40 capemaxcels 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 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) 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 SAVE nid_day, nid_mth, nid_ins 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 REAL t2m(klon,nbsrf), q2m(klon,nbsrf) REAL u10m(klon,nbsrf), v10m(klon,nbsrf) REAL zt2m(klon), zq2m(klon) REAL zu10m(klon), zv10m(klon) CHARACTER*40 t2mincels, t2maxcels 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 histISCCP #include "ini_histISCCP.h" #endif #ifdef histhf #include "ini_histhf.h" #endif #include "ini_histday.h" #include "ini_histmth.h" #include "ini_histins.h" 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 c IF(rmu0(i).EQ.0.) THEN c sunlit(i)=0 c PRINT*,' il fait nuit ',i,rlat(i),rlon(i) c ENDIF 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, 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 ISCCP simulator BEGIN IF (ok_isccp) THEN cIM calcul tau. emi nuages convectifs convfra(:,:)=rnebcon(:,:) convliq(:,:)=rnebcon(:,:)*clwcon(:,:) c CALL newmicro (paprs, pplay,ok_newmicro, c . t_seri, cldliq, cldfra, cldtau, cldemi, c . cldh, cldl, cldm, cldt, cldq) CALL newmicro (paprs, pplay,ok_newmicro, . t_seri, convliq, convfra, dtau_c, dem_c, . cldh_c, cldl_c, cldm_c, cldt_c, cldq_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) cIM calcul diagramme (PC, tau) cf. ISCCP D c seed=50 c seed=ran0(klon) cT1O3 c top_height=1 cT3O3 c top_height=3 c overlap=3 cIM cf GCM cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.) cIM inversion des niveaux de pression ==> de haut en bas DO k=1,klev kinv=klev-k+1 DO i=1,klon pfull(i,k)=pplay(i,kinv) c on met toutes les variables de Haut 2 Bas qv(i,k)=q_seri(i,kinv) cc(i,k)=cldtot(i,kinv) conv(i,k)=rnebcon(i,kinv) dtau_sH2B(i,k)=dtau_s(i,kinv) dtau_cH2B(i,k)=dtau_c(i,kinv) at(i,k)=t_seri(i,kinv) dem_sH2B(i,k)=dem_s(i,kinv) dem_cH2B(i,k)=dem_c(i,kinv) ENDDO ENDDO DO k=1,klev+1 kinv=klev-k+2 DO i=1,klon phalf(i,k)=paprs(i,kinv) ENDDO ENDDO 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 nsrf=3 DO nreg=1, nbreg 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 c pct_ocean(i,nreg)=.FALSE. pct_ocean(i,nreg)=0 c DO nsrf = 1, nbsrf 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 c pct_ocean(i,nreg)=.TRUE. 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 c pct_ocean(i,nreg)=.TRUE. 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 c pct_ocean(i,nreg)=.TRUE. 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 c pct_ocean(i,nreg)=.TRUE. 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 c pct_ocean(i,nreg)=.TRUE. pct_ocean(i,nreg)=1 ENDIF ENDIF ENDIF !nbreg 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 c ENDDO ENDDO !klon ENDDO !nbreg cIM somme de toutes les nhistoW BEG DO nreg = 1, nbreg DO k = 1, kmaxm1 DO l = 1, lmaxm1 DO iw = 1, iwmax nhistoWt(k,l,iw,nreg)=0. ENDDO ENDDO ENDDO ENDDO cIM somme de toutes les nhistoW END ENDIF c CALL ISCCP_CLOUD_TYPES(nlev,ncol,seed,pfull,phalf,qv, c & cc,conv,dtau_s,dtau_c,top_height,overlap, c & tautab,invtau,skt,emsfc_lw,at,dem_s,dem_c,fq_isccp, c & totalcldarea,meanptop,meantaucld,boxtau,boxptop) c DO i=1, klon c i=1 c1011 CONTINUE c cIM on verifie les donnees de INPUT en dehors du simulateur ISCCP cIM 1D non-vectorise (!) pour qu'on gagne du temps ... cIM c BEGIN find unpermittable data..... ! ---------------------------------------------------! ! find unpermittable data..... ! c do 13 k=1,klev c ca prend trop de temps ?? c cldtot(:,:) = min(max(cldtot(:,:),0.),1.) c rnebcon(:,:) = min(max(rnebcon(:,:),0.),1.) c dtau_s(:,:) = max(dtau_s(:,:),0.) c dem_s(:,:) = min(max(dem_s(:,:),0.),1.) c dtau_c(:,:) = max(dtau_c(:,:),0.) c dem_c(:,:) = min(max(dem_c(:,:),0.),1.) c ca prend trop de temps ?? c if (cldtot(i,k) .lt. 0.) then c print *, ' error = cloud fraction less than zero' c STOP c end if c if (cldtot(i,k) .gt. 1.) then c print *, ' error = cloud fraction greater than 1' c STOP c end if c if (rnebcon(i,k) .lt. 0.) then c print *, c & ' error = convective cloud fraction less than zero' c STOP c end if c if (rnebcon(i,k) .gt. 1.) then c print *, c & ' error = convective cloud fraction greater than 1' c STOP c end if c if (dtau_s(i,k) .lt. 0.) then c print *, c & ' error = stratiform cloud opt. depth less than zero' c STOP c end if c if (dem_s(i,k) .lt. 0.) then c print *, c & ' error = stratiform cloud emissivity less than zero' c STOP c end if c if (dem_s(i,k) .gt. 1.) then c print *, c & ' error = stratiform cloud emissivity greater than 1' c STOP c end if c if (dtau_c(i,k) .lt. 0.) then c print *, c & ' error = convective cloud opt. depth less than zero' c STOP c end if c if (dem_c(i,k) .lt. 0.) then c print *, c & ' error = convective cloud emissivity less than zero' c STOP c end if c if (dem_c(i,k) .gt. 1.) then c print *, c & ' error = convective cloud emissivity greater than 1' c STOP c end if c13 continue ! ---------------------------------------------------! c c END find unpermittable data..... cv2.2.1.1 DO i=1, klon c i=1 c seed=i c cv3.4 if (debut) then DO i=1, klon seed(i)=i+100 c seed(i)=i+50 ENDDO endif c seed=aint(ran0(klon)) c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed,pfull(i,:),phalf(i,:) cv2.2.1.1 c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed(i),pfull(i,:),phalf(i,:) c & ,q_seri(i,:), c & cldtot(i,:),rnebcon(i,:),dtau_s(i,:),dtau_c(i,:), c & top_height,overlap, c & tautab,invtau,ztsol,emsfc_lw,t_seri(i,:),dem_s(i,:), c & dem_c(i,:), c & fq_isccp(i,:,:), c & totalcldarea(i),meanptop(i),meantaucld(i), c & boxtau(i,:),boxptop(i,:)) cv2.2.1.1 cv3.4 debug=0 debugcol=0 cIM260503 c o500 ==> distribution nuage ftion du regime dynamique DO i=1, klon o500(i)=omega(i,8)*864. c PRINT*,'pphi8 ',pphi(i,8),'zphi8,11,12',zphi(i,8), c & zphi(i,11),zphi(i,12) ENDDO c axe vertical pour les differents niveaux des histogrammes c DO iw=1, iwmax c zx_o500(iw)=wmin+(iw-1./2.)*pas_w c ENDDO c PRINT*,' phys AVANT seed(3361)=',seed(3361) CALL ISCCP_CLOUD_TYPES( & debug, & debugcol, & klon, & sunlit, & klev, & ncol, & seed, & pfull, & phalf, c var de bas en haut ==> PB ! c & q_seri, c & cldtot, c & rnebcon, c & dtau_s, c & dtau_c, c var de Haut en Bas BEG & qv, cc, conv, dtau_sH2B, dtau_cH2B, c var de Haut en Bas END & top_height, & overlap, & tautab, & invtau, & ztsol, & emsfc_lw, c var de bas en haut ==> PB ! c & t_seri, c & dem_s, c & dem_c, c var de Haut en Bas BEG & at, dem_sH2B, dem_cH2B, cIM260503 c & o500, pct_ocean, c var de Haut en Bas END & fq_isccp, & totalcldarea, & meanptop, & meantaucld, & boxtau, & boxptop) c & boxptop, cIM 260503 c & histoW, c & nhistoW c &) cIM 200603 c PRINT*,'physiq fq_isccp(6,1,1)',fq_isccp(6,1,1) cIM 200603 cIM somme de toutes les nhistoW BEG c DO k = 1, kmaxm1 c DO l = 1, lmaxm1 c DO iw = 1, iwmax c nhistoWt(k,l,iw)=nhistoWt(k,l,iw)+nhistoW(k,l,iw) ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then c IF(nhistoWt(k,l,iw).NE.0.) THEN c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw) c ENDIF c ENDDO c ENDDO c ENDDO cIM somme de toutes les nhistoW END c PRINT*,' phys APRES seed(3361)=',seed(3361) cv3.4 c i=i+1 c IF(i.LE.klon) THEN c GOTO 1011 c ENDIF cv2.2.1.1 ENDDO c passage de la grille (klon,7,7) a (iim,jjmp1,7,7) c minfq3d=100. c maxfq3d=0. cIM calcul des correspondances entre la grille physique et cIM la grille dynamique c DO i=1, klon c grid_phys(i)=i c PRINT*,'i, grid_phys',i,grid_phys(i) c ENDDO c CALL gr_fi_dyn(1,klon,iimp1,jjmp1,grid_phys,grid_dyn) c DO j=1, jjmp1 c DO i=1, iimp1 c PRINT*,'i,j grid_dyn ',i,j,grid_dyn(i,j) c ENDDO c ENDDO c DO l=1, lmax DO k=1, kmax cIM grille physique ==> grille ecriture 2D (iim,jjmp1) c DO i=1, iim fq4d(i,1,k,l)=fq_isccp(1,k,l) cc PRINT*,'first j=1 i =',i ENDDO DO j=2, jjm DO i=1, iim cERROR ?? ig=i+iim*(j-1) ig=i+1+(j-2)*iim cc PRINT*,'i =',i,'j =',j,'ig=',ig 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) cc PRINT*,'last jjmp1 i =',i ENDDO IF(debut) THEN DO j=1, jjmp1 DO i=1, iim IF(j.GE.2.AND.j.LE.jjm) THEN igfi2D(i,j)=i+1+(j-2)*iim c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j) ELSEIF(j.EQ.1) THEN igfi2D(i,j)=1 c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j) ELSEIF(j.EQ.jjmp1) THEN igfi2D(i,j)=klon c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j) ENDIF ENDDO ENDDO ENDIF c STOP c c CALL gr_fi_ecrit(1,klon,iim,jjmp1,fq_isccp(:,k,l), c $ fq4d(:,:,k,l)) ENDDO ENDDO DO l=1, lmax fq4d(:,:,8,l)=-1.e+10 fq4d(:,:,l,8)=-1.e+10 ENDDO DO l=1, lmax DO k=1, kmax DO j=1, jjmp1 DO i=1, iim c inversion TAU ?! c ni=(i-1)*lmax+l c nj=(j-1)*kmax+kmax-k+1 c c210503 inversion en PC ==> pas besoin !!! c ni=(i-1)*lmax+lmax-l+1 c nj=(j-1)*kmax+k c c210503 ni=(i-1)*lmax+l nj=(j-1)*kmax+k c210503 if(k.EQ.8) then c fq4d(i,j,8,l)=-1.e+10 c endif c210503 if(l.EQ.8) then c fq4d(i,j,k,8)=-1.e+10 c endif fq3d(ni,nj)=fq4d(i,j,k,l) c if(fq3d(ni,nj).lt.0.) then c print*,' fq3d LT ZERO ',ni,nj,fq3d(ni,nj) c endif c if(fq3d(ni,nj).gt.100.) then c print*,' fq3d GT 100 ',ni,nj,fq3d(ni,nj) c endif c max & min fq3d c if(fq3d(ni,nj).gt.maxfq3d) maxfq3d=fq3d(ni,nj) c if(fq3d(ni,nj).lt.minfq3d) minfq3d=fq3d(ni,nj) ENDDO ENDDO c fq4d(:,:,8,l)=-1.e+10 c fq4d(:,:,k,8)=-1.e+10 c k=k+1 c if(k.LE.kmax) then c goto 1022 c endif ENDDO c l=l+1 c if(l.LE.lmax) then c goto 1021 c endif ENDDO c print*,' minfq3d=',minfq3d,' maxfq3d=',maxfq3d c c calculs statistiques distribution nuage ftion du regime dynamique c DO i=1, klon c! o500(i)=omega(i,9)*864. c! PRINT*,' o500=',o500(i),' pphi(9)=',pphi(i,9) c o500(i)=omega(i,8)*864. cc PRINT*,' pphi(8)',pphi(i,8),'pphi(11)',pphi(i,11), cc .'pphi(12)',pphi(i,12) cc PRINT*,' zphi8,11,12=',zphi(i,8),zphi(i,11),zphi(i,12) cc PRINT*,' o500',o500(i),' w500',w500(i) c ENDDO c axe vertical pour les differents niveaux des histogrammes c DO iw=1, iwmax c zx_o500(iw)=wmin+(iw-1./2.)*pas_w c ENDDO c Ce calcul doit etre fait a partir de valeurs mensuelles ?? cc CALL histo_o500_pctau(o500,fq4d,histoW) cc CALL histo_o500_pctau(paire,pctsrf,o500,fq4d,histoW) cc CALL histo_o500_pctau(pct_ocean,rlat,o500,fq4d,histoW) ccOK ??? CALL histo_o500_pctau(pct_ocean,o500,fq4d,histoW) c CALL histo_o500_pctau(klon,pct_ocean,o500,fq4d,histoW,nhistoW) c CALL histo_o500_pctau(klon,pct_ocean,o500,fq_isccp, CALL histo_o500_pctau(nbreg,pct_ocean,o500,fq_isccp, &histoW,nhistoW) c cIM somme de toutes les nhistoW BEG DO nreg=1, nbreg 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) ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then c IF(nhistoWt(k,l,iw).NE.0.) THEN c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw) c ENDIF ENDDO ENDDO ENDDO ENDDO cIM somme de toutes les nhistoW END c c IF(lafin) THEN c DO nreg=1, nbreg c DO iw=1, iwmax c DO l=1,lmaxm1 c DO k=1,kmaxm1 c IF(histoW(k,l,iw,nreg).NE.0.) then c PRINT*,'physiq H nH',k,l,iw, c & histoW(k,l,iw,nreg), c & nhistoW(k,l,iw,nreg),nhistoWt(k,l,iw,nreg) c ENDIF c ENDDO c ENDDO c ENDDO c ENDDO cIM verif fq_isccp, fq4d, fq3d c DO l=1, lmaxm1 c DO k=1,kmaxm1 c i=74 c j=36 c DO j=1, jjmp1 c DO i=1, iim c DO l=1, lmaxm1 c WRITE(*,'(a,3i4,7f10.4)') c & 'fq_isccp,j,i,l=',j,i,l, c & (fq_isccp(igfi2D(i,j),k,l),k=1,kmaxm1) c WRITE(*,'(a,3i4,7f10.4)') c & 'fq4d,j,i,l=',j,i,l,(fq4d(i,j,k,l),k=1,kmaxm1) c ENDDO c ENDDO c ENDDO c ni1=(i-1)*8+1 c ni2=i*8 c nj1=(j-1)*8+1 c nj2=j*8 c DO ni=ni1,ni2 c WRITE(*,'(a,2i4,7f10.4)') c & 'fq3d, ni,nj=',ni,nj, c & (fq3d(ni,nj),nj=nj1,nj2) c ENDDO c ENDIF c DO iw=1, iwmax c DO l=1,lmaxm1 c DO k=1,kmaxm1 c PRINT*,' iw,l,k,nhistoW=',iw,l,k,nhistoW(k,l,iw) c ENDDO c ENDDO c ENDDO c DO iw=1, iwmax c DO l=1, lmaxm1 c linv=lmaxm1-l+1 c DO k=1, kmaxm1 c histoWinv(k,l,iw)=histoW(iw,k,l) c ENDDO c ENDDO c ENDDO c c pb syncronisation ?? : 48 * 30 * 7 (jour1) + 48* 29 * 7 (jour suivant) c ENDIF !ok_isccp cIM ISCCP simulator END 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) 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) cIM e (dist, rmu0, fract, co2_ppm, solaire, 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, cccIMs topsw0,toplw0,solsw0,sollw0) 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 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 cccIM cf. FH c======================================================================= c SORTIES c======================================================================= c Interpollation sur quelques niveaux de pression c ----------------------------------------------- call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850) call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850) call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500) call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500) call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200) call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200) call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500) call plevel(klon,klev,.true. ,paprs,50000.,omega,w500) cIM cf. FH slp(:) = paprs(:,1)*exp(pphis(:)/(289.*t_seri(:,1))) slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1))) c PRINT*,' physiq slp ',slp(2185),paprs(2185,1),pphis(2185), c . RD,t_seri(2185,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 c PRINT*,' i ',i,' prw',prw(i) ENDDO c c============================================================= c Ecriture des sorties c============================================================= #ifdef histISCCP #include "write_histISCCP.h" #endif #ifdef histhf #include "write_histhf.h" #endif #include "write_histday.h" #include "write_histmth.h" #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