Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F (revision 503) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/physiq.F (revision 504) @@ -124,7 +124,8 @@ c PARAMETER (ok_journe=.true.) - integer lev_histday - save lev_histday - data lev_histday/1/ +cIM lev_histday ==> clesphys +cIM integer lev_histday +cIM save lev_histday +cIM data lev_histday/1/ c LOGICAL ok_mensuel ! sortir le fichier mensuel @@ -195,26 +196,57 @@ #include "raddim.h" c - REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2) +cIM 080304 REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2) + REAL swdn0(klon,klevp1), swdn(klon,klevp1) + REAL swup0(klon,klevp1), swup(klon,klevp1) SAVE swdn0 , swdn, swup0, swup - +c + REAL SWdn200clr(klon), SWdn200(klon) + REAL SWup200clr(klon), SWup200(klon) + SAVE SWdn200clr, SWdn200, SWup200clr, SWup200 +c + REAL lwdn0(klon,klevp1), lwdn(klon,klevp1) + REAL lwup0(klon,klevp1), lwup(klon,klevp1) + SAVE lwdn0 , lwdn, lwup0, lwup +c + REAL LWdn200clr(klon), LWdn200(klon) + REAL LWup200clr(klon), LWup200(klon) + SAVE LWdn200clr, LWdn200, LWup200clr, LWup200 +c + REAL LWdnTOA(klon), LWdnTOAclr(klon) + SAVE LWdnTOA, LWdnTOAclr +c 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 + integer nlevSTD + PARAMETER(nlevSTD=17) + real rlevSTD(nlevSTD) + DATA rlevSTD/100000., 92500., 85000., 70000., + .60000., 50000., 40000., 30000., 25000., 20000., + .15000., 10000., 7000., 5000., 3000., 2000., 1000./ + CHARACTER*5 clevSTD(nlevSTD), aa, bb + DATA clevSTD/'1000','925 ','850 ','700 ','600 ', + .'500 ','400 ','300 ','250 ','200 ','150 ','100 ', + .'70 ','50 ','30 ','20 ','10 '/ +c + real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD) + real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD) + real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD) +c +cIM ENSEMBLES BEG +c + integer nlevENS + PARAMETER(nlevENS=4) + integer indENS(nlevENS) + save indENS + real rlevENS(nlevENS) + DATA rlevENS/85000., 70000., 50000., 20000./ + CHARACTER*3 clev(nlevENS) + DATA clev/'850','700','500','200'/ + + real tlev(klon,nlevENS), qlev(klon,nlevENS), rhlev(klon,nlevENS) + real ulev(klon,nlevENS), vlev(klon,nlevENS), philev(klon,nlevENS) + real wlev(klon,nlevENS) +cIM ENSEMBLES END +c c prw: precipitable water real prw(klon) @@ -322,4 +354,5 @@ SAVE nhistoWt + INTEGER linv INTEGER pct_ocean(klon,nbregdyn) REAL rlonPOS(klon) @@ -380,9 +413,10 @@ logical ok_hf real ecrit_hf - integer nid_hf - save ok_hf, ecrit_hf, nid_hf + integer nid_hf, nid_hf3d + save ok_hf, ecrit_hf, nid_hf, nid_hf3d c QUESTION : noms de variables ? +#undef histhf #define histhf #ifdef histhf @@ -566,4 +600,8 @@ REAL snow_fall(klon) ! neige save snow_fall, rain_fall +cIM 050204 BEG + REAL total_rain(klon), nday_rain(klon) + save total_rain, nday_rain +cIM 050204 END REAL evap(klon), devap(klon) ! evaporation et sa derivee REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee @@ -571,5 +609,5 @@ 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 +C ! type de sous-surface et pondere par la fraction REAL fder(klon) ! Derive de flux (sensible et latente) save fder @@ -590,4 +628,8 @@ SAVE lmt_pas ! frequence de mise a jour REAL pctsrf(klon,nbsrf) +cIM + REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE + REAL paire_ter(klon) !surfaces terre +cIM SAVE pctsrf ! sous-fraction du sol REAL albsol(klon) @@ -632,4 +674,7 @@ cIM EXTERNAL haut2bas !variables de haut en bas + INTEGER lnblnk1 + EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type + !caracter c c Variables locales @@ -664,7 +709,11 @@ REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon) real sollwdown(klon) ! downward LW flux at surface +cIM BEG + real sollwdownclr(klon) ! downward CS LW flux at surface + real toplwdown(klon) ! downward CS LW flux at TOA + real toplwdownclr(klon) ! downward CS LW flux at TOA +cIM END 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 @@ -672,7 +721,7 @@ c sauvegarder les sorties du rayonnement SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown + SAVE sollwdownclr, toplwdown, toplwdownclr SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0 -cccIM - +c INTEGER itaprad SAVE itaprad @@ -832,5 +881,5 @@ c INTEGER nhori, nvert - REAL zsto, zout + REAL zsto, zout, zsto1, zsto2 real zjulian save zjulian @@ -899,8 +948,34 @@ IF (debut) THEN CALL suphec ! initialiser constantes et parametres phys. - ENDIF - - +c +cIM 050204 BEG + DO i=1, klon + nday_rain(i)=0. + ENDDO +cIM 050204 END +c c====================================================================== +cIM BEG + DO k=1, nlev + DO l=1, nlevSTD +c + bb=clevSTD(l) +c + IF(l.GE.2) THEN + aa=clevSTD(l) + bb=aa(1:lnblnk1(aa)) + ENDIF +c + IF(bb.EQ.clev(k)) THEN +c print*,'k=',k,'l=',l,'clev=',clev(k) + indENS(k)=l +c print*,'k=',k,'l=',l,'clev=',clev(k),'indENS=',indENS(k) + ENDIF +c + ENDDO + ENDDO +c + ENDIF !debut +cIM END xjour = rjourvrai c @@ -927,5 +1002,4 @@ itap = 0 itaprad = 0 -c CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0, . rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow, @@ -938,5 +1012,5 @@ radpas = NINT( 86400./dtime/nbapp_rad) c -C on remet le calendrier à zero +C on remet le calendrier a zero c IF (raz_date .eq. 1) THEN @@ -1053,5 +1127,5 @@ c c -cccIM +cIM capemaxcels = 't_max(X)' t2mincels = 't_min(X)' @@ -1156,12 +1230,12 @@ ENDIF C + DO i = 1, klon + ztsol(i) = 0. + ENDDO + DO nsrf = 1, nbsrf DO i = 1, klon - ztsol(i) = 0. + ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf) ENDDO - DO nsrf = 1, nbsrf - DO i = 1, klon - ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf) - ENDDO - ENDDO + ENDDO C IF (if_ebil.ge.1) THEN @@ -1324,5 +1398,7 @@ fder = dlw - CALL clmain(dtime,itap,date0,pctsrf, + +cIM CALL clmain(dtime,itap,date0,pctsrf, + CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new, e t_seri,q_seri,u_seri,v_seri, e julien, rmu0, co2_ppm, @@ -1339,5 +1415,4 @@ 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) @@ -1399,5 +1474,5 @@ zxtsol(i) = 0.0 zxfluxlat(i) = 0.0 -cccIM +c zt2m(i) = 0.0 zq2m(i) = 0.0 @@ -1952,10 +2027,10 @@ ENDDO !nreg cIM somme de toutes les nhistoW END -c + ENDIF cIM: initialisation de seed DO i=1, klon seed(i)=i+100 ENDDO - ENDIF !debut + cIM: pas de debug, debugcol debug=0 @@ -2152,5 +2227,10 @@ s topsw,toplw,solsw,sollw, s sollwdown, +cIM s sollwdown, sollwdownclr, +cIM +cIM s toplwdown, toplwdownclr, +cIM s topsw0,toplw0,solsw0,sollw0, + s lwdn0, lwdn, lwup0, lwup, s swdn0, swdn, swup0, swup ) itaprad = 0 @@ -2393,59 +2473,101 @@ 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 on moyenne mensuellement les champs 3D et on les interpole sur les niveaux STD +c if(itap.EQ.1.OR.itap.EQ.13.OR.itap.EQ.25.OR.itap.EQ.37) THEN +c if(MOD(itap,12).EQ.1) THEN +cIM 120304 END + DO k=1, nlevSTD + call plevel(klon,klev,.true.,pplay,rlevSTD(k), + . t_seri,tlevSTD(:,k)) + call plevel(klon,klev,.false.,pplay,rlevSTD(k), + . u_seri,ulevSTD(:,k)) + call plevel(klon,klev,.false.,pplay,rlevSTD(k), + . v_seri,vlevSTD(:,k)) + call plevel(klon,klev,.false.,pplay,rlevSTD(k), + . zphi,philevSTD(:,k)) + call plevel(klon,klev,.false.,pplay,rlevSTD(k), + . qx(:,:,ivap),qlevSTD(:,k)) + call plevel(klon,klev,.false.,pplay,rlevSTD(k), + . zx_rh,rhlevSTD(:,k)) + ENDDO !nlevSTD +c ENSEMBLES BEG + DO k=1, nlevENS +cIM 170304 + tlev(:,k)=tlevSTD(:,indENS(k)) + ulev(:,k)=ulevSTD(:,indENS(k)) + vlev(:,k)=vlevSTD(:,indENS(k)) + philev(:,k)=philevSTD(:,indENS(k)) + qlev(:,k)=qlevSTD(:,indENS(k)) + rhlev(:,k)=rhlevSTD(:,indENS(k)) +c + call plevel(klon,klevp1,.true.,paprs,rlevENS(k), + . omega,wlev(:,k)) +c + ENDDO !k=1, nlevENS +cIM 170304 +c IF(1.EQ.0) THEN +c +c call plevel(klon,klev,.true.,pplay,rlevENS(k), +c . t_seri,tlev(:,k)) +c . t_seri,tlevSTD(:,indENS(k))) +c call plevel(klon,klev,.false.,pplay,rlevENS(k), +c . qx(:,:,ivap),qlev(:,k)) +c . qx(:,:,ivap),qlevSTD(:,indENS(k))) +c call plevel(klon,klev,.false.,pplay,rlevENS(k), +c . zx_rh,rhlev(:,k)) +c . zx_rh,rhlevSTD(:,indENS(k))) +c +c IF(k.EQ.1) THEN +c ulev(:,k)=u850(:) +c vlev(:,k)=v850(:) +c philev(:,k)=phi850(:) +c ulev(:,k)=ulevSTD(:,3) +c vlev(:,k)=vlevSTD(:,3) +c philev(:,k)=philevSTD(:,3) +c ELSEIF(k.EQ.2) THEN +c ulev(:,k)=u500(:) +c vlev(:,k)=v500(:) +c philev(:,k)=phi500(:) +c ulev(:,k)=ulevSTD(:,6) +c vlev(:,k)=vlevSTD(:,6) +c philev(:,k)=philevSTD(:,6) +c ELSEIF(k.EQ.3) THEN +c ulev(:,k)=u200(:) +c vlev(:,k)=v200(:) +c philev(:,k)=phi200(:) +c ulev(:,k)=ulevSTD(:,10) +c vlev(:,k)=vlevSTD(:,10) +c philev(:,k)=philevSTD(:,10) +c ENDIF !k +c +c ENDIF !1.EQ.0 +c +c ENDDO !nlevENS +c120304 ENDIF !IF (MOD(itap,ecrit_mth) .EQ. 0) THEN +c endif !(itap.EQ.ecrit_mth) THEN +cIM 100304 BEG +cIM interpolation a chaque pas de temps du SWup(clr) et SWdn(clr) a 200 hPa + call plevel(klon,klevp1,.true.,paprs,20000., + $ swdn0,SWdn200clr) + call plevel(klon,klevp1,.false.,paprs,20000., + $ swdn,SWdn200) + call plevel(klon,klevp1,.false.,paprs,20000., + $ swup0,SWup200clr) + call plevel(klon,klevp1,.false.,paprs,20000., + $ swup,SWup200) +c + call plevel(klon,klevp1,.false.,paprs,20000., + $ lwdn0,LWdn200clr) + call plevel(klon,klevp1,.false.,paprs,20000., + $ lwdn,LWdn200) + call plevel(klon,klevp1,.false.,paprs,20000., + $ lwup0,LWup200clr) + call plevel(klon,klevp1,.false.,paprs,20000., + $ lwup,LWup200) +c +cIM 100304 END +c +c ENSEMBLES END +c c slp sea level pressure slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1))) @@ -2467,4 +2589,59 @@ ENDDO ENDDO +cIM 050204 BEG +c + IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN +cIM PRINT *,' PHYS cond julien ',julien +c CALL ozonecm( FLOAT(julien), rlat, paprs, wo) + DO i = 1, klon + total_rain(i)=rain_fall(i)+snow_fall(i) + IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1. + ENDDO +c + ENDIF +c surface terre + IF (debut) THEN + DO i=1, klon + IF(pctsrf_new(i,is_ter).GT.0.) THEN + paire_ter(i)=paire(i)*pctsrf_new(i,is_ter) + ENDIF + ENDDO + ENDIF +cIM 050204 END + +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 22.03.04 BEG c============================================================= c Ecriture des sorties @@ -2485,42 +2662,11 @@ #include "write_histmth.h" -#ifdef histmthNMC -#include "write_histmthNMC.h" -#endif +c#ifdef histmthNMC +c#include "write_histmthNMC.h" +c#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 22.03.04 END c c====================================================================