Index: trunk/LMDZ.MARS/libf/dyn3d/iniprint.h =================================================================== --- trunk/LMDZ.MARS/libf/dyn3d/iniprint.h (revision 1535) +++ (revision ) @@ -1,11 +1,0 @@ -! -! $Id: $ -! -! -! handle debug and output levels -! lunout: unit of file where outputs will be sent -! (default is 6, standard output) -! prt_level: level of informative output messages (0 = minimum) -! - INTEGER lunout, prt_level - COMMON /comprint/ lunout, prt_level Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/calfis.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/calfis.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/calfis.F (revision 1540) @@ -0,0 +1,554 @@ + SUBROUTINE calfis(nq, lafin, rdayvrai,rday_ecri, heure, + $ pucov,pvcov,pteta,pq,pmasse,pps,pp,ppk,pphis,pphi, + $ pducov,pdvcov,pdteta,pdq,pw, + $ pdufi,pdvfi,pdhfi,pdqfi,pdpsfi,tracer ) +c +c Auteur : P. Le Van, F. Hourdin +c ......... + + USE comvert_mod, ONLY: preff + USE comconst_mod, ONLY: dtphys,cpp,kappa,pi + + IMPLICIT NONE +c======================================================================= +c +c 1. rearrangement des tableaux et transformation +c variables dynamiques > variables physiques +c 2. calcul des termes physiques +c 3. retransformation des tendances physiques en tendances dynamiques +c +c remarques: +c ---------- +c +c - les vents sont donnes dans la physique par leurs composantes +c naturelles. +c - la variable thermodynamique de la physique est une variable +c intensive : T +c pour la dynamique on prend T * ( preff / p(l) ) **kappa +c - les deux seules variables dependant de la geometrie necessaires +c pour la physique sont la latitude pour le rayonnement et +c l'aire de la maille quand on veut integrer une grandeur +c horizontalement. +c - les points de la physique sont les points scalaires de la +c la dynamique; numerotation: +c 1 pour le pole nord +c (jjm-1)*iim pour l'interieur du domaine +c ngridmx pour le pole sud +c ---> ngridmx=2+(jjm-1)*iim +c +c Input : +c ------- +c ecritphy frequence d'ecriture (en jours)de histphy +c pucov covariant zonal velocity +c pvcov covariant meridional velocity +c pteta potential temperature +c pps surface pressure +c pmasse masse d'air dans chaque maille +c pts surface temperature (K) +c pw flux vertical (kg/s) +c +c Output : +c -------- +c pdufi tendency for the natural zonal velocity (ms-1) +c pdvfi tendency for the natural meridional velocity +c pdhfi tendency for the potential temperature +c pdtsfi tendency for the surface temperature +c +c tracer Call tracer in gcm.F ? (decided in callphys.def) +c +c======================================================================= +c +c----------------------------------------------------------------------- +c +c 0. Declarations : +c ------------------ + +#include "dimensions.h" +#include "paramet.h" + + INTEGER ngridmx,nq + PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) + +#include "comgeom2.h" +!#include "control.h" + +c Arguments : +c ----------- + LOGICAL lafin + REAL heure + + REAL pvcov(iip1,jjm,llm) + REAL pucov(iip1,jjp1,llm) + REAL pteta(iip1,jjp1,llm) + REAL pmasse(iip1,jjp1,llm) + REAL pq(iip1,jjp1,llm,nq) + REAL pphis(iip1,jjp1) + REAL pphi(iip1,jjp1,llm) +c + REAL pdvcov(iip1,jjm,llm) + REAL pducov(iip1,jjp1,llm) + REAL pdteta(iip1,jjp1,llm) + REAL pdq(iip1,jjp1,llm,nq) +c + REAL pw(iip1,jjp1,llm) +c + REAL pps(iip1,jjp1) + REAL pp(iip1,jjp1,llmp1) + REAL ppk(iip1,jjp1,llm) +c + REAL pdvfi(iip1,jjm,llm) + REAL pdufi(iip1,jjp1,llm) + REAL pdhfi(iip1,jjp1,llm) + REAL pdqfi(iip1,jjp1,llm,nq) + REAL pdpsfi(iip1,jjp1) + logical tracer + +c Local variables : +c ----------------- + + INTEGER i,j,l,ig0,ig,iq + REAL zpsrf(ngridmx) + REAL zplev(ngridmx,llm+1),zplay(ngridmx,llm) + REAL zphi(ngridmx,llm),zphis(ngridmx) +c + REAL zufi(ngridmx,llm), zvfi(ngridmx,llm) + REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nq) +c + REAL zvervel(ngridmx,llm) +c + REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm) + REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nq) + REAL zdpsrf(ngridmx) +c + REAL zsin(iim),zcos(iim),z1(iim) + REAL zsinbis(iim),zcosbis(iim),z1bis(iim) + REAL unskap, pksurcp +c + + EXTERNAL gr_dyn_fi,gr_fi_dyn + EXTERNAL physiq,multipl + REAL SSUM + EXTERNAL SSUM + + REAL latfi(ngridmx),lonfi(ngridmx) + REAL airefi(ngridmx) + SAVE latfi, lonfi, airefi + + LOGICAL firstcal, debut + DATA firstcal/.true./ + SAVE firstcal,debut + REAL rdayvrai,rday_ecri +c +c----------------------------------------------------------------------- +c +c 1. Initialisations : +c -------------------- +c + + + IF (ngridmx.NE.2+(jjm-1)*iim) THEN + PRINT*,'STOP dans calfis' + PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim' + PRINT*,' ngridmx jjm iim ' + PRINT*,ngridmx,jjm,iim + STOP + ENDIF + +c----------------------------------------------------------------------- +c latitude, longitude et aires des mailles pour la physique: +c ---------------------------------------------------------- + +c + IF ( firstcal ) THEN + debut = .TRUE. + ELSE + debut = .FALSE. + ENDIF + +c +! IF (firstcal) THEN +! latfi(1)=rlatu(1) +! lonfi(1)=0. +! DO j=2,jjm +! DO i=1,iim +! latfi((j-2)*iim+1+i)= rlatu(j) +! lonfi((j-2)*iim+1+i)= rlonv(i) +! ENDDO +! ENDDO +! latfi(ngridmx)= rlatu(jjp1) +! lonfi(ngridmx)= 0. +! +! ! build airefi(), mesh area on physics grid +! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) +! ! Poles are single points on physics grid +! airefi(1)=airefi(1)*iim +! airefi(ngridmx)=airefi(ngridmx)*iim +! +! CALL inifis(ngridmx,llm,nq,day_ini,daysec,dtphys, +! . latfi,lonfi,airefi,rad,g,r,cpp) +! ENDIF + +c +c----------------------------------------------------------------------- +c 40. transformation des variables dynamiques en variables physiques: +c --------------------------------------------------------------- + +c 41. pressions au sol (en Pascals) +c ---------------------------------- + + + zpsrf(1) = pps(1,1) + + ig0 = 2 + DO j = 2,jjm + CALL SCOPY( iim,pps(1,j),1,zpsrf(ig0), 1 ) + ig0 = ig0+iim + ENDDO + + zpsrf(ngridmx) = pps(1,jjp1) + + +c 42. pression intercouches : +c +c ----------------------------------------------------------------- +c .... zplev definis aux (llm +1) interfaces des couches .... +c .... zplay definis aux ( llm ) milieux des couches .... +c ----------------------------------------------------------------- + +c ... Exner = cp * ( p(l) / preff ) ** kappa .... +c + unskap = 1./ kappa +c + DO l = 1, llmp1 + zplev( 1,l ) = pp(1,1,l) + ig0 = 2 + DO j = 2, jjm + DO i =1, iim + zplev( ig0,l ) = pp(i,j,l) + ig0 = ig0 +1 + ENDDO + ENDDO + zplev( ngridmx,l ) = pp(1,jjp1,l) + ENDDO +c +c + +c 43. temperature naturelle (en K) et pressions milieux couches . +c --------------------------------------------------------------- + + DO l=1,llm + + pksurcp = ppk(1,1,l) / cpp + zplay(1,l) = preff * pksurcp ** unskap + ztfi(1,l) = pteta(1,1,l) * pksurcp + ig0 = 2 + + DO j = 2, jjm + DO i = 1, iim + pksurcp = ppk(i,j,l) / cpp + zplay(ig0,l) = preff * pksurcp ** unskap + ztfi(ig0,l) = pteta(i,j,l) * pksurcp + ig0 = ig0 + 1 + ENDDO + ENDDO + + pksurcp = ppk(1,jjp1,l) / cpp + zplay(ig0,l) = preff * pksurcp ** unskap + ztfi (ig0,l) = pteta(1,jjp1,l) * pksurcp + + ENDDO + + DO l=1, llm + DO ig=1,ngridmx + if (ztfi(ig,l).lt.15) then + write(*,*) 'New Temperature below 15 K !!! ' + write(*,*) 'Stop in calfis.F ' + write(*,*) 'ig=', ig, ' l=', l + write(*,*) 'ztfi(ig,l)=',ztfi(ig,l) + stop + end if + ENDDO + ENDDO + + + +c 43.bis Taceurs (en kg/kg) +c -------------------------- + DO iq=1,nq + DO l=1,llm + zqfi(1,l,iq) = pq(1,1,l,iq) + ig0 = 2 + DO j=2,jjm + DO i = 1, iim + zqfi(ig0,l,iq) = pq(i,j,l,iq) + ig0 = ig0 + 1 + ENDDO + ENDDO + zqfi(ig0,l,iq) = pq(1,jjp1,l,iq) + ENDDO + ENDDO + +c Geopotentiel calcule par rapport a la surface locale: +c ----------------------------------------------------- + + CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,pphi,zphi) + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,pphis,zphis) + DO l=1,llm + DO ig=1,ngridmx + zphi(ig,l)=zphi(ig,l)-zphis(ig) + ENDDO + ENDDO + +c Calcul de la vitesse verticale (m/s) pour diagnostique +c ------------------------------- +c pw est en kg/s +c On interpole "lineairement" la temperature entre les couches(FF,10/95) + +! DO ig=1,ngridmx +! zvervel(ig,1)=0. +! END DO +! DO l=2,llm +! zvervel(1,l)=(pw(1,1,l)/apoln) +! & * r *0.5*(ztfi(1,l)+ztfi(1,l-1)) /zplev(1,l) +! ig0=2 +! DO j=2,jjm +! DO i = 1, iim +! zvervel(ig0,l) = pw(i,j,l) * unsaire(i,j) +! & * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l) +! ig0 = ig0 + 1 +! ENDDO +! ENDDO +! zvervel(ig0,l)=(pw(1,jjp1,l)/apols) +! & * r *0.5*(ztfi(ig0,l)+ztfi(ig0,l-1)) /zplev(ig0,l) +! ENDDO + +c ......... Reindexation : calcul de zvervel au MILIEU des couches +! DO l=1,llm-1 +! DO ig=1,ngridmx +! zvervel(ig,l) = 0.5*(zvervel(ig,l)+zvervel(ig,l+1)) +! END DO +! END DO +c (dans la couche llm, on garde la valeur à la limite inférieure llm) + +c 45. champ u: +c ------------ + + DO 50 l=1,llm + + DO 25 j=2,jjm + ig0 = 1+(j-2)*iim + zufi(ig0+1,l)= 0.5 * + $ ( pucov(iim,j,l)/cu(iim,j) + pucov(1,j,l)/cu(1,j) ) + DO 10 i=2,iim + zufi(ig0+i,l)= 0.5 * + $ ( pucov(i-1,j,l)/cu(i-1,j) + pucov(i,j,l)/cu(i,j) ) +10 CONTINUE +25 CONTINUE + +50 CONTINUE + + +c 46.champ v: +c ----------- + + DO l=1,llm + DO j=2,jjm + ig0=1+(j-2)*iim + DO i=1,iim + zvfi(ig0+i,l)= 0.5 * + $ ( pvcov(i,j-1,l)/cv(i,j-1) + pvcov(i,j,l)/cv(i,j) ) + ENDDO + ENDDO + ENDDO + + +c 47. champs de vents aux pole nord +c ------------------------------ +c U = 1 / pi * integrale [ v * cos(long) * d long ] +c V = 1 / pi * integrale [ v * sin(long) * d long ] + + DO l=1,llm + + z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,1,l)/cv(1,1) + z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,1,l)/cv(1,1) + DO i=2,iim + z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,1,l)/cv(i,1) + z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,1,l)/cv(i,1) + ENDDO + + DO i=1,iim + zcos(i) = COS(rlonv(i))*z1(i) + zcosbis(i)= COS(rlonv(i))*z1bis(i) + zsin(i) = SIN(rlonv(i))*z1(i) + zsinbis(i)= SIN(rlonv(i))*z1bis(i) + ENDDO + + zufi(1,l) = SSUM(iim,zcos,1)/pi + zvfi(1,l) = SSUM(iim,zsin,1)/pi + + ENDDO + + +c 48. champs de vents aux pole sud: +c --------------------------------- +c U = 1 / pi * integrale [ v * cos(long) * d long ] +c V = 1 / pi * integrale [ v * sin(long) * d long ] + + DO l=1,llm + + z1(1) =(rlonu(1)-rlonu(iim)+2.*pi)*pvcov(1,jjm,l)/cv(1,jjm) + z1bis(1)=(rlonu(1)-rlonu(iim)+2.*pi)*pdvcov(1,jjm,l)/cv(1,jjm) + DO i=2,iim + z1(i) =(rlonu(i)-rlonu(i-1))*pvcov(i,jjm,l)/cv(i,jjm) + z1bis(i)=(rlonu(i)-rlonu(i-1))*pdvcov(i,jjm,l)/cv(i,jjm) + ENDDO + + DO i=1,iim + zcos(i) = COS(rlonv(i))*z1(i) + zcosbis(i) = COS(rlonv(i))*z1bis(i) + zsin(i) = SIN(rlonv(i))*z1(i) + zsinbis(i) = SIN(rlonv(i))*z1bis(i) + ENDDO + + zufi(ngridmx,l) = SSUM(iim,zcos,1)/pi + zvfi(ngridmx,l) = SSUM(iim,zsin,1)/pi + + ENDDO + +c----------------------------------------------------------------------- +c Appel de la physique: +c --------------------- + + + CALL physiq (ngridmx,llm,nq, + , debut,lafin, + , rday_ecri,heure,dtphys, + , zplev,zplay,zphi, + , zufi, zvfi,ztfi, zqfi, +! , zvervel, + , pw, +C - sorties + s zdufi, zdvfi, zdtfi, zdqfi,zdpsrf,tracer) + + +c----------------------------------------------------------------------- +c transformation des tendances physiques en tendances dynamiques: +c --------------------------------------------------------------- + +c tendance sur la pression : +c ----------------------------------- + + CALL gr_fi_dyn(1,ngridmx,iip1,jjp1,zdpsrf,pdpsfi) +c +ccc CALL multipl(ip1jmp1,aire,pdpsfi,pdpsfi) + +c 62. enthalpie potentielle +c --------------------- + + DO l=1,llm + + DO i=1,iip1 + pdhfi(i,1,l) = cpp * zdtfi(1,l) / ppk(i, 1 ,l) + pdhfi(i,jjp1,l) = cpp * zdtfi(ngridmx,l)/ ppk(i,jjp1,l) + ENDDO + + DO j=2,jjm + ig0=1+(j-2)*iim + DO i=1,iim + pdhfi(i,j,l) = cpp * zdtfi(ig0+i,l) / ppk(i,j,l) + ENDDO + pdhfi(iip1,j,l) = pdhfi(1,j,l) + ENDDO + + ENDDO + + +c 62. humidite specifique +c --------------------- + + DO iq=1,nq + DO l=1,llm + DO i=1,iip1 + pdqfi(i,1,l,iq) = zdqfi(1,l,iq) + pdqfi(i,jjp1,l,iq) = zdqfi(ngridmx,l,iq) + ENDDO + DO j=2,jjm + ig0=1+(j-2)*iim + DO i=1,iim + pdqfi(i,j,l,iq) = zdqfi(ig0+i,l,iq) + ENDDO + pdqfi(iip1,j,l,iq) = pdqfi(1,j,l,iq) + ENDDO + ENDDO + ENDDO + +c 65. champ u: +c ------------ + + DO l=1,llm + + DO i=1,iip1 + pdufi(i,1,l) = 0. + pdufi(i,jjp1,l) = 0. + ENDDO + + DO j=2,jjm + ig0=1+(j-2)*iim + DO i=1,iim-1 + pdufi(i,j,l)= + $ 0.5*(zdufi(ig0+i,l)+zdufi(ig0+i+1,l))*cu(i,j) + ENDDO + pdufi(iim,j,l)= + $ 0.5*(zdufi(ig0+1,l)+zdufi(ig0+iim,l))*cu(iim,j) + pdufi(iip1,j,l)=pdufi(1,j,l) + ENDDO + + ENDDO + + +c 67. champ v: +c ------------ + + DO l=1,llm + + DO j=2,jjm-1 + ig0=1+(j-2)*iim + DO i=1,iim + pdvfi(i,j,l)= + $ 0.5*(zdvfi(ig0+i,l)+zdvfi(ig0+i+iim,l))*cv(i,j) + ENDDO + pdvfi(iip1,j,l) = pdvfi(1,j,l) + ENDDO + ENDDO + + +c 68. champ v pres des poles: +c --------------------------- +c v = U * cos(long) + V * SIN(long) + + DO l=1,llm + + DO i=1,iim + pdvfi(i,1,l)= + $ zdufi(1,l)*COS(rlonv(i))+zdvfi(1,l)*SIN(rlonv(i)) + pdvfi(i,jjm,l)=zdufi(ngridmx,l)*COS(rlonv(i)) + $ +zdvfi(ngridmx,l)*SIN(rlonv(i)) + pdvfi(i,1,l)= + $ 0.5*(pdvfi(i,1,l)+zdvfi(i+1,l))*cv(i,1) + pdvfi(i,jjm,l)= + $ 0.5*(pdvfi(i,jjm,l)+zdvfi(ngridmx-iip1+i,l))*cv(i,jjm) + ENDDO + + pdvfi(iip1,1,l) = pdvfi(1,1,l) + pdvfi(iip1,jjm,l)= pdvfi(1,jjm,l) + + ENDDO + +c----------------------------------------------------------------------- + +700 CONTINUE + + firstcal = .FALSE. + + RETURN + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_dyn_fi.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_dyn_fi.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_dyn_fi.F (revision 1540) @@ -0,0 +1,37 @@ + SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi) + + IMPLICIT NONE +c======================================================================= +c passage d'un champ de la grille scalaire a la grille physique +c======================================================================= + +c----------------------------------------------------------------------- +c declarations: +c ------------- + + INTEGER im,jm,ngrid,nfield + REAL pdyn(im,jm,nfield) + REAL pfi(ngrid,nfield) + + INTEGER j,ifield,ig + EXTERNAL SCOPY + +c----------------------------------------------------------------------- +c calcul: +c ------- + + IF(ngrid.NE.2+(jm-2)*(im-1)) STOP 'probleme de dim' +c traitement des poles + CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid) + CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid) + +c traitement des point normaux + DO ifield=1,nfield + DO j=2,jm-1 + ig=2+(j-2)*(im-1) + CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1) + ENDDO + ENDDO + + RETURN + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_fi_dyn.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_fi_dyn.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/gr_fi_dyn.F (revision 1540) @@ -0,0 +1,38 @@ + SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn) + IMPLICIT NONE +c======================================================================= +c passage d'un champ de la grille scalaire a la grille physique +c======================================================================= + +c----------------------------------------------------------------------- +c declarations: +c ------------- + + INTEGER im,jm,ngrid,nfield + REAL pdyn(im,jm,nfield) + REAL pfi(ngrid,nfield) + + INTEGER i,j,ifield,ig + EXTERNAL SCOPY + +c----------------------------------------------------------------------- +c calcul: +c ------- + + DO ifield=1,nfield +c traitement des poles + DO i=1,im + pdyn(i,1,ifield)=pfi(1,ifield) + pdyn(i,jm,ifield)=pfi(ngrid,ifield) + ENDDO + +c traitement des point normaux + DO j=2,jm-1 + ig=2+(j-2)*(im-1) + CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1) + pdyn(im,j,ifield)=pdyn(1,j,ifield) + ENDDO + ENDDO + + RETURN + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/caldyn0.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/caldyn0.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/caldyn0.F (revision 1540) @@ -0,0 +1,1 @@ +link ../../dyn3d/caldyn0.F Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F (revision 1540) @@ -0,0 +1,301 @@ +c======================================================================= + SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0, + & zmea,zstd,zsig,zgam,zthe) +c======================================================================= +c +c +c Author: F. Hourdin 01/1997 +c ------- +c +c Object: To read data from Martian surface to use in a GCM +c ------ from NetCDF file "surface.nc" +c +c +c Arguments: +c ---------- +c +c Inputs: +c ------ +c +c Outputs: +c -------- +c +c======================================================================= +c donnees ALBEDO, INERTIE THERMIQUE, RELIEF: +c +c Ces donnees sont au format NetCDF dans le fichier "surface.nc" +c +c 360 valeurs en longitude (de -179.5 a 179.5) +c 180 valeurs en latitudes (de 89.5 a -89.5) +c +c Pour les passer au format de la grille, on utilise "interp_horiz.F" +c +c Il faut donc que ces donnees soient au format grille scalaire +c (imold+1 jmold+1) +c avec passage des coordonnees de la "boite" (rlonu, rlatv) +c +c On prend imd (d pour donnees!) +c imd = 360 avec copie de la 1ere valeur sur la imd+1 +c (rlonud de -179 a -181) +c jmd = 179 +c (rlatvd de 89 a -89) +c======================================================================= + +! to use 'getin' + use ioipsl_getincom + USE comconst_mod, ONLY: g,pi + + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom.h" +#include "netcdf.inc" +#include "datafile.h" + +c======================================================================= +c Declarations: +C======================================================================= + + INTEGER imd,jmd,imdp1,jmdp1 + parameter (imd=360,jmd=179,imdp1=361,jmdp1=180) + + INTEGER iimp1 + parameter (iimp1=iim+1-1/iim) + +! Arguments: + CHARACTER(len=3),intent(inout) :: relief + REAL,intent(out) :: phisinit(iimp1*jjp1) + REAL,intent(out) :: alb(iimp1*jjp1) + REAL,intent(out) :: ith(iimp1*jjp1) + REAL,intent(out) :: z0(iimp1*jjp1) + REAL,intent(out) :: zmea(imdp1*jmdp1) + REAL,intent(out) :: zstd(imdp1*jmdp1) + REAL,intent(out) :: zsig(imdp1*jmdp1) + REAL,intent(out) :: zgam(imdp1*jmdp1) + REAL,intent(out) :: zthe(imdp1*jmdp1) + +! Local variables: + REAL zdata(imd*jmdp1) + REAL zdataS(imdp1*jmdp1) + REAL pfield(iimp1*jjp1) + + INTEGER ierr + + INTEGER unit,nvarid + + INTEGER i,j,k + + INTEGER klatdat,ngridmxgdat + PARAMETER (klatdat=180,ngridmxgdat=360) + +c on passe une grille en rlonu rlatv et im+1 jm a interp_horiz) + + REAL longitude(imd),latitude(jmdp1) ! Pour lecture des donnees + REAL rlonud(imdp1),rlatvd(jmd) + + CHARACTER*20 string + DIMENSION string(0:4) + + +!#include "lmdstd.h" +!#include "fxyprim.h" + + pi=2.*ASIN(1.) + +c======================================================================= +c rlonud, rlatvd +c======================================================================= + +c----------------------------------------------------------------------- +c Lecture NetCDF des donnees latitude et longitude +c----------------------------------------------------------------------- + write(*,*) 'datareadnc: opening file surface.nc' + + datafile="/u/lmdz/WWW/planets/mars/datadir" ! default path to surface.nc + call getin("datadir",datafile) ! but users may specify another path + + ierr = NF_OPEN (trim(datafile)//'/surface.nc', + & NF_NOWRITE,unit) + IF (ierr.NE.NF_NOERR) THEN + write(*,*)'Error : cannot open file surface.nc ' + write(*,*)'(in phymars/datareadnc.F)' + write(*,*)'It should be in :',trim(datafile),'/' + write(*,*)'1) You can set this path in the + & callphys.def file:' + write(*,*)' datadir=/path/to/the/datafiles' + write(*,*)'2) If necessary, surface.nc (and other datafiles)' + write(*,*)' can be obtained online on:' + write(*,*)'http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' + CALL ABORT + ENDIF + +c +c Lecture des latitudes (coordonnees): +c + ierr = NF_INQ_VARID (unit, "latitude", nvarid) +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(unit, nvarid, latitude) +#else + ierr = NF_GET_VAR_REAL(unit, nvarid, latitude) +#endif +c +c Lecture des longitudes (coordonnees): +c + ierr = NF_INQ_VARID (unit, "longitude", nvarid) +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(unit, nvarid, longitude) +#else + ierr = NF_GET_VAR_REAL(unit, nvarid, longitude) +#endif + +c----------------------------------------------------------------------- +c Passage au format boites scalaires +c----------------------------------------------------------------------- + +c----------------------------------------------------------------------- +c longitude(imd) --> rlonud(imdp1) +c----------------------------------------------------------------------- + +c Passage en coordonnees boites scalaires et en radian + do i=1,imd + rlonud(i)=(longitude(i)+.5)*pi/180. + enddo + +c Repetition de la valeur im+1 + rlonud(imdp1)=rlonud(1) + 2*pi + +c----------------------------------------------------------------------- +c latitude(jmdp1) --> rlonvd(jmd) +c----------------------------------------------------------------------- + +c Passage en coordonnees boites scalaires et en radian + do j=1,jmd + rlatvd(j)=(latitude(j)-.5)*pi/180. + enddo + +c======================================================================= +c lecture NetCDF de albedo, thermal, relief, zdtm (pour francois Lott) +c======================================================================= + + string(0) = 'z0' + string(1) = 'albedo' + string(2) = 'thermal' + if (relief.ne.'pla') then + write(*,*) ' MOLA topography' + relief = 'MOL' + string(3) = 'z'//relief + else + string(3) = 'zMOL' ! pour qu''il lise qqchose sur le fichier + ! remise a 0 derriere + endif + string(4) = 'zMOL' ! lecture pour calcul topog. sous-maille + + + DO k=0,4 + write(*,*) 'string',k,string(k) + +c----------------------------------------------------------------------- +c initialisation +c----------------------------------------------------------------------- + call initial0(iimp1*jjp1,pfield) + call initial0(imd*jmdp1,zdata) + call initial0(imdp1*jmdp1,zdataS) + +c----------------------------------------------------------------------- +c Lecture NetCDF +c----------------------------------------------------------------------- + + ierr = NF_INQ_VARID (unit, string(k), nvarid) + if (ierr.ne.nf_noerr) then + write(*,*) 'datareadnc error, cannot find ',trim(string(k)) + write(*,*) ' in file ',trim(datafile),'/surface.nc' + stop + endif +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata) +#else + ierr = NF_GET_VAR_REAL(unit, nvarid, zdata) +#endif + if (ierr.ne.nf_noerr) then + write(*,*) 'datareadnc: error failed loading ',trim(string(k)) + stop + endif + +c----------------------------------------------------------------------- +c Cas particulier "Francois Lott" ( k=4 ) (relief sous-maille) +c----------------------------------------------------------------------- + if (k.eq.4) then + + zdata(:)=1000.*zdata(:) + longitude(:)=(pi/180.)*longitude(:) + latitude(:)=(pi/180.)*latitude(:) + + call grid_noro1(360, 180, longitude, latitude, zdata, + . iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe) + + endif + +c----------------------------------------------------------------------- +c Passage de zdata en grille (imdp1 jmdp1) +c----------------------------------------------------------------------- + do j=1,jmdp1 + do i=1,imd + zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1)) + enddo + zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1)) + enddo + +c----------------------------------------------------------------------- +c Interpolation +c----------------------------------------------------------------------- + call interp_horiz(zdataS,pfield,imd,jmd, + . iim, jjm,1,rlonud,rlatvd,rlonu,rlatv) + +c----------------------------------------------------------------------- +c Periodicite +c----------------------------------------------------------------------- + + do j=1,jjp1 + pfield(iimp1*j) = pfield(1+iimp1*(j-1)) + enddo + +c----------------------------------------------------------------------- +c Sauvegarde des champs +c----------------------------------------------------------------------- + + if (k.eq.0) then ! z0 + z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01 + ! multiplied by 0.01 to have z0 in m + elseif (k.eq.1) then ! albedo + do i=1,iimp1*jjp1 + alb(i) = pfield(i) + enddo + elseif (k.eq.2) then ! thermal + do i=1,iimp1*jjp1 + ith(i) = pfield(i) + enddo + elseif (k.eq.3) then ! relief + if (relief.eq.'pla') then + call initial0(iimp1*jjp1,phisinit) + else + do i=1,iimp1*jjp1 + phisinit(i) = pfield(i) + enddo + endif + endif + + ENDDO + +c----------------------------------------------------------------------- +c Traitement Phisinit +c----------------------------------------------------------------------- + + phisinit(1:iimp1*jjp1)=1000.*phisinit(1:iimp1*jjp1) + phisinit(:)=g*phisinit(:) + +c----------------------------------------------------------------------- +c FIN +c----------------------------------------------------------------------- + + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/defrun_new.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/defrun_new.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/defrun_new.F (revision 1540) @@ -0,0 +1,1 @@ +link ../../dyn3d/defrun_new.F Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/exner_hyb.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/exner_hyb.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/exner_hyb.F (revision 1540) @@ -0,0 +1,1 @@ +link ../../dyn3d/exner_hyb.F Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/grid_noro1.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/grid_noro1.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/grid_noro1.F (revision 1540) @@ -0,0 +1,425 @@ + SUBROUTINE grid_noro1(imdep, jmdep, xdata, ydata, entree, + . imar, jmar, x, y, zmea,zstd,zsig,zgam,zthe) +c======================================================================= +c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead) +c +c Calcul des parametres de l'orographie sous-maille necessaires +c au nouveau shema de representation des montagnes meso-echelles +c dans le modele. Les points sont mis sur une grille rectangulaire +c pseudo-physique. Typiquement, il y a iim+1 latitudes incluant +c le pole nord et le pole sud. Il y a jjm+1 longitudes, y compris +c aux poles. Aux poles les champs peuvent ont une valeurs repetee +c jjm+1 fois..... La valeur du champs en jjm+1 (jmar) est celle +c en j=1. +c Les parametres a,b,c,d representent les limites de la region +c de point de grille correspondant a un point decrit precedemment. +c Les moyennes sur ces regions des valeurs calculees a partir de +c l'USN, sont ponderees par un poids, fonction de la surface +c occuppe par ces donnees a l'interieure de la grille du modele. +c Dans la plupart des cas ce poid est le rapport entre la surface +c de la region de point de grille USN et la surface de la region +c de point de grille du modele. +c +c +c (c) +c ----d----- +c | . . . .| +c | | +c (b)a . * . .b(a) +c | | +c | . . . .| +c ----c----- +c (d) +C======================================================================= +c INPUT: +c imdep, jmdep: dimensions X et Y pour depart +c xdata, ydata: coordonnees X et Y pour depart +c entree: champ d'entree a transformer +c dans ce programme, on assume que les donnees sont les altitudes +c de l'USNavy: imdep=iusn=2160, jmdep=jusn=1080. +c OUTPUT: +c imar, jmar: dimensions X et Y d'arrivee +c x, y: coordonnees X et Y d'arrivee +c les champs de sorties sont sur une grille physique: +c zmea: orographie moyenne +c zstd: deviation standard de l'orographie sous-maille +c zsig: pente de l'orographie sous-maille +c zgam: anisotropy de l'orographie sous maille +c zthe: orientation de l'axe oriente dans la direction +c de plus grande pente de l'orographie sous maille +C======================================================================= +c IMPLICIT INTEGER (I,J) +c IMPLICIT REAL(X,Z) + + USE comconst_mod, ONLY: rad + + implicit none + integer iusn,jusn,iext + parameter(iusn=360,jusn=180,iext=40) +c!-*- include 'param1' +c!-*- include 'comcstfi.h' +#include "dimensions.h" +c!-*- +c!-*- parameter(iim=cols,jjm=rows) + REAL xusn(iusn+2*iext),yusn(jusn+2) + REAL zusn(iusn+2*iext,jusn+2),zusnfi(iusn+2*iext,jusn+2) + +c modif declarations pour implicit none + real zmeanor,zmeasud,zstdnor,zstdsud,zsignor + real zsigsud,zweinor,zweisud + real xk,xl,xm,xw,xp,xq + real zmaxmea,zmaxstd,zmaxsig,zmaxgam,zmaxthe,zminthe + real zbordnor,zbordsud,zbordest,zbordoue,xpi + real zdeltax,zdeltay,zlenx,zleny,weighx,weighy,xincr + integer i,j,ii,jj,ideltax,ihalph + + INTEGER imdep, jmdep + REAL xdata(imdep),ydata(jmdep) + REAL entree(imdep,jmdep) +c + INTEGER imar, jmar + + REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1) + REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1) + REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2) + REAL zxtzyusn(iusn+2*iext,jusn+2) + REAL weight(iim+1,jjm+1) + REAL x(imar+1),y(jmar) + REAL zmea(imar+1,jmar),zstd(imar+1,jmar) + REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar) +c + REAL a(2200),b(2200),c(1100),d(1100) +c +c quelques constantes: +c + print *,' parametres de l orographie a l echelle sous maille' + print*,'rad =',rad + print*,'Long et lat entree' + print*,(x(i),i=1,imar+1) + print*,(y(j),j=1,jmar) + print*,'Long et lat donnees' + print*,(xdata(i),i=1,imdep) + print*,(ydata(j),j=1,jmdep) + + xpi=acos(-1.) + zdeltay=2.*xpi/real(jusn)*rad +c +c quelques tests de dimensions: +c + IF (imar.GT.2200 .OR. jmar.GT.1100) THEN + PRINT*, 'imar ou jmar trop grand', imar, jmar + CALL ABORT + ENDIF + + IF(imdep.ne.iusn.or.jmdep.ne.jusn)then + print *,' imdep ou jmdep mal dimensionnes:',imdep,jmdep + call abort + ENDIF + + IF(imar+1.gt.iim+1.or.jmar.gt.jjm+1)THEN + print *,' imar ou jmar mal dimensionnes:',imar,jmar + call abort + ENDIF +c +C Extension de la base de donnee de l'USN pour faciliter +C les calculs ulterieurs: +c + DO j=1,jusn + yusn(j+1)=ydata(j) + DO i=1,iusn + zusn(i+iext,j+1)=entree(i,j) + xusn(i+iext)=xdata(i) + ENDDO + DO i=1,iext + zusn(i,j+1)=entree(iusn-iext+i,j) + xusn(i)=xdata(iusn-iext+i)-2.*xpi + zusn(iusn+iext+i,j+1)=entree(i,j) + xusn(iusn+iext+i)=xdata(i)+2.*xpi + ENDDO + ENDDO + + yusn(1)=ydata(1)+(ydata(1)-ydata(2)) + yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1)) + DO i=1,iusn/2+iext + zusn(i,1)=zusn(i+iusn/2,2) + zusn(i+iusn/2+iext,1)=zusn(i,2) + zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1) + zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1) + ENDDO +c +c Calcul d'une orographie filtree aux hautes latitudes +c pour permettre des calculs plus isotropiques sur la pente +c des montagnes +c + DO i=1,IUSN+2*iext + DO J=1,JUSN+2 + zusnfi(i,j)=0.0 + ENDDO + ENDDO + + DO j=1,jusn + ideltax=1./cos(yusn(j+1)) + ideltax=min(iusn/2-1,ideltax) + IF(MOD(IDELTAX,2).EQ.0)THEN + IDELTAX=IDELTAX+1 + ENDIF + IHALPH=(IDELTAX-1)/2 +c print *,' ideltax=',ideltax + IF(ideltax.eq.1)THEN + DO i=1,iusn + zusnfi(i+iext,j+1)=entree(i,j) + ENDDO + ELSE + DO i=1,ihalph + DO ii=1,i+ihalph + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j) + ENDDO + DO ii=ihalph-i,0,-1 + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(iusn-ii,j) + ENDDO + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax) + ENDDO + DO i=iusn-ihalph+1,iusn + DO ii = i-ihalph,iusn + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j) + ENDDO + DO ii = 1,ihalph+i-iusn + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(ii,j) + ENDDO + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax) + ENDDO + DO i=ihalph+1,iusn-ihalph + DO ii=-ihalph,ihalph + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)+entree(i+ii,j) + ENDDO + zusnfi(i+iext,j+1)=zusnfi(i+iext,j+1)/real(ideltax) + ENDDO + ENDIF + DO i=1,iext + zusnfi(i,j+1)=zusnfi(iusn-iext+i,j+1) + zusnfi(i+iusn+iext,j+1)=zusnfi(i,j+1) + ENDDO + ENDDO +c +c Calculer les limites des zones des nouveaux points +c + a(1) = x(1) - (x(2)-x(1))/2.0 + b(1) = (x(1)+x(2))/2.0 + DO i = 2, imar-1 + a(i) = b(i-1) + b(i) = (x(i)+x(i+1))/2.0 + ENDDO + a(imar) = b(imar-1) + b(imar) = x(imar) + (x(imar)-x(imar-1))/2.0 + + c(1) = y(1) - (y(2)-y(1))/2.0 + d(1) = (y(1)+y(2))/2.0 + DO j = 2, jmar-1 + c(j) = d(j-1) + d(j) = (y(j)+y(j+1))/2.0 + ENDDO + c(jmar) = d(jmar-1) + d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0 +c +c quelques initialisations: + print*,'OKM1' +c + DO i = 1, imar + DO j = 1, jmar + weight(i,j) = 0.0 + zxtzx(i,j) = 0.0 + zytzy(i,j) = 0.0 + zxtzy(i,j) = 0.0 + ztz(i,j) = 0.0 + zmea(i,j) = 0.0 + zstd(i,j)=0.0 + ENDDO + ENDDO +c +c calculs des correlations de pentes sur la grille de l'USN. +c + DO j = 2,jusn+1 + DO i = 1, iusn+2*iext + zytzyusn(i,j)=0.0 + zxtzxusn(i,j)=0.0 + zxtzyusn(i,j)=0.0 + ENDDO + ENDDO + + + DO j = 2,jusn+1 + zdeltax=zdeltay*cos(yusn(j)) + DO i = 2, iusn+2*iext-1 + zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2 + zxtzxusn(i,j)=(zusnfi(i+1,j)-zusnfi(i-1,j))**2/zdeltax**2 + zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay + * *(zusnfi(i+1,j)-zusnfi(i-1,j))/zdeltax + ENDDO + + ENDDO + + + + print*,'OK0' +c +c sommations des differentes quantites definies precedemment +c sur une grille du modele. +c + zleny=xpi/real(jusn)*rad + xincr=xpi/2./real(jusn) + DO ii = 1, imar + DO jj = 1, jmar +c PRINT *,' iteration ii jj:',ii,jj + DO j = 2,jusn+1 +c DO j = 3,jusn + zlenx=zleny*cos(yusn(j)) + zdeltax=zdeltay*cos(yusn(j)) + zbordnor=(c(jj)-yusn(j)+xincr)*rad + zbordsud=(yusn(j)-d(jj)+xincr)*rad + weighy=amax1(0., + * amin1(zbordnor,zbordsud,zleny)) + IF(weighy.ne.0)THEN + DO i = 2, iusn+2*iext-1 + zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j)) + zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j)) + weighx=amax1(0., + * amin1(zbordest,zbordoue,zlenx)) + IF(weighx.ne.0)THEN + weight(ii,jj)=weight(ii,jj)+weighx*weighy + zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy + zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy + zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy + ztz(ii,jj) =ztz(ii,jj) +zusn(i,j)*zusn(i,j)*weighx*weighy + zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy + ENDIF + ENDDO + ENDIF + ENDDO + ENDDO + ENDDO +c +c calculs des differents parametres necessaires au programme +c de parametrisation de l'orographie a l'echelle moyenne: +c + zmaxmea=0. + zmaxstd=0. + zmaxsig=0. + zmaxgam=0. + zmaxthe=0. + zminthe=0. +c print 100,' ' +c100 format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH') + print*,'OK1' + DO ii = 1, imar + DO jj = 1, jmar +c print*,'ok0' + IF (weight(ii,jj) .NE. 0.0) THEN +c Orography moyenne: +c print*,'ok1' + zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj) + zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj) + zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj) + zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj) + ztz(ii,jj) =ztz(ii,jj)/weight(ii,jj) +c print*,'ok2' +c Deviation standard: + zstd(ii,jj)=sqrt(amax1(0.,ztz(ii,jj)-zmea(ii,jj)**2)) +c Coefficients K, L et M: + xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2. + xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2. + xm=zxtzy(ii,jj) + xp=xk-sqrt(xl**2+xm**2) + xq=xk+sqrt(xl**2+xm**2) + xw=1.e-8 + if(xp.le.xw) xp=0. + if(xq.le.xw) xq=xw + if(abs(xm).le.xw) xm=xw*sign(1.,xm) +c print*,'ok3' +c pente: + zsig(ii,jj)=sqrt(xq) +c zsig(ii,jj)=sqrt(2.*xk) +c isotropy: + zgam(ii,jj)=xp/xq +c angle theta: + zthe(ii,jj)=57.29577951*atan2(xm,xl)/2. + +c print 101,ii,jj, +c * zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj), +c * zthe(ii,jj) +c101 format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1) +c print*,'ok4' + ELSE +c PRINT*, 'probleme,ii,jj=', ii,jj +c print*,'ok1b' + ENDIF + zmaxmea=amax1(zmea(ii,jj),zmaxmea) +c print*,'oka' + zmaxstd=amax1(zstd(ii,jj),zmaxstd) +c print*,'okb' + zmaxsig=amax1(zsig(ii,jj),zmaxsig) +c print*,'okc' + zmaxgam=amax1(zgam(ii,jj),zmaxgam) +c print*,'okd' + zmaxthe=amax1(zthe(ii,jj),zmaxthe) +c print*,'oke' + zminthe=amin1(zthe(ii,jj),zminthe) +c print*,'ok5' + ENDDO + ENDDO + + print *,' MEAN ORO:',zmaxmea + print *,' ST. DEV.:',zmaxstd + print *,' PENTE:',zmaxsig + print *,' ANISOTROP:',zmaxgam + print *,' ANGLE:',zminthe,zmaxthe + +C +c On passe ce donnees sur la grille dite physique....(?) +c On met gamma et theta a 1. et 0. aux poles ou ces quantites +c n'ont pas vraiment de sens +c + DO jj=1,jmar + zmea(imar+1,jj)=zmea(1,jj) + zstd(imar+1,jj)=zstd(1,jj) + zsig(imar+1,jj)=zsig(1,jj) + zgam(imar+1,jj)=zgam(1,jj) + zthe(imar+1,jj)=zthe(1,jj) + ENDDO + + + zmeanor=0.0 + zmeasud=0.0 + zstdnor=0.0 + zstdsud=0.0 + zsignor=0.0 + zsigsud=0.0 + zweinor=0.0 + zweisud=0.0 + + DO ii=1,imar + zweinor=zweinor+ weight(ii, 1) + zweisud=zweisud+ weight(ii,jmar) + zmeanor=zmeanor+zmea(ii, 1)*weight(ii, 1) + zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar) + zstdnor=zstdnor+zstd(ii, 1)*weight(ii, 1) + zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar) + zsignor=zsignor+zsig(ii, 1)*weight(ii, 1) + zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar) + ENDDO + + DO ii=1,imar+1 + zmea(ii, 1)=zmeanor/zweinor + zmea(ii,jmar)=zmeasud/zweisud + zstd(ii, 1)=zstdnor/zweinor + zstd(ii,jmar)=zstdsud/zweisud + zsig(ii, 1)=zsignor/zweinor + zsig(ii,jmar)=zsigsud/zweisud + zgam(ii, 1)=1. + zgam(ii,jmar)=1. + zthe(ii, 1)=0. + zthe(ii,jmar)=0. + ENDDO + + + RETURN + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/ini_archive.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/ini_archive.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/ini_archive.F (revision 1540) @@ -0,0 +1,521 @@ +c======================================================================= + subroutine ini_archive(nid,idayref,phis,ith,tab_cntrl_fi) +c======================================================================= +c +c +c Date: 01/1997 +c ---- +c +c Objet: ecriture de l'entete du fichier "start_archive" +c ----- +c +c Proche de iniwrite.F +c +c On ajoute dans le tableau "tab_cntrl" (dynamique), a partir de 51, +c les valeurs de tab_cntrl_fi (les 38 parametres de controle physiques +c du RUN + ptotal et cotoicetotal) +c +c tab_cntrl(50+l)=tab_cntrl_fi(l) +c +c Arguments: +c --------- +c +c Inputs: +c ------ +c +c nid unite logique du fichier "start_archive" +c idayref Valeur du jour initial a mettre dans +c l'entete du fichier "start_archive" +c phis geopotentiel au sol +c ith soil thermal inertia +c tab_cntrl_fi tableau des param physiques +c + +c======================================================================= + + use comsoil_h, only: nsoilmx, mlayer + USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,presnivs,pseudoalt + USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi + USE logic_mod, ONLY: fxyhypb,ysinus + USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy + USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom.h" +#include "netcdf.inc" + +c----------------------------------------------------------------------- +c Declarations +c----------------------------------------------------------------------- + +c Local: +c ------ + INTEGER length,l + parameter (length = 100) + REAL tab_cntrl(length) ! tableau des parametres du run + INTEGER loop + INTEGER ierr, setvdim, putvdim, putdat, setname,cluvdb + INTEGER setdim + INTEGER ind1,indlast + +c Arguments: +c ---------- + INTEGER*4 idayref + REAL phis(ip1jmp1) + real ith(ip1jmp1,nsoilmx) + REAL tab_cntrl_fi(length) + +!Mars --------Ajouts----------- +c Variables locales pour NetCDF: +c + INTEGER dims2(2), dims3(3) !, dims4(4) + INTEGER idim_index + INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_rlatv + INTEGER idim_llmp1,idim_llm + INTEGER idim_tim + INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID # + INTEGER nid,nvarid + real sig_s(llm),s(llm) + + pi = 2. * ASIN(1.) + + +c----------------------------------------------------------------------- +c Remplissage du tableau des parametres de controle du RUN (dynamique) +c----------------------------------------------------------------------- + + DO l=1,length + tab_cntrl(l)=0. + ENDDO + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + tab_cntrl(1) = REAL(iim) ! nombre de points en longitude + tab_cntrl(2) = REAL(jjm) ! nombre de points en latitude + tab_cntrl(3) = REAL(llm) ! nombre de couches + tab_cntrl(4) = REAL(idayref) ! jour 0 + tab_cntrl(5) = rad ! rayon de mars(m) ~3397200 + tab_cntrl(6) = omeg ! vitesse de rotation (rad.s-1) + tab_cntrl(7) = g ! gravite (m.s-2) ~3.72 + tab_cntrl(8) = cpp + tab_cntrl(8) = 43.49 !mars temporaire Masse molaire de l''atm (g.mol-1) ~43.49 + tab_cntrl(9) = kappa ! = r/cp ~0.256793 (=rcp dans physique) + tab_cntrl(10) = daysec ! duree du sol (s) ~88775 + tab_cntrl(11) = dtvr ! pas de temps de la dynamique (s) + tab_cntrl(12) = etot0 ! energie totale ! + tab_cntrl(13) = ptot0 ! pression totalei ! variables + tab_cntrl(14) = ztot0 ! enstrophie totale ! de controle + tab_cntrl(15) = stot0 ! enthalpie totale ! globales + tab_cntrl(16) = ang0 ! moment cinetique ! + tab_cntrl(17) = pa + tab_cntrl(18) = preff + +c ..... parametres pour le zoom ...... + + tab_cntrl(19) = clon ! longitude en degres du centre du zoom + tab_cntrl(20) = clat ! latitude en degres du centre du zoom + tab_cntrl(21) = grossismx ! facteur de grossissement du zoom,selon longitude + tab_cntrl(22) = grossismy ! facteur de grossissement du zoom ,selon latitude + + IF ( fxyhypb ) THEN + tab_cntrl(23) = 1. + tab_cntrl(24) = dzoomx ! extension en longitude de la zone du zoom + tab_cntrl(25) = dzoomy ! extension en latitude de la zone du zoom + ELSE + tab_cntrl(23) = 0. + tab_cntrl(24) = dzoomx ! extension en longitude de la zone du zoom + tab_cntrl(25) = dzoomy ! extension en latitude de la zone du zoom + tab_cntrl(26) = 0. + IF ( ysinus) tab_cntrl(26) = 1. + ENDIF + +c----------------------------------------------------------------------- +c Copie du tableau des parametres de controle du RUN (physique) +c dans le tableau dynamique +c----------------------------------------------------------------------- + + DO l=1,50 + tab_cntrl(50+l)=tab_cntrl_fi(l) + ENDDO + +c======================================================================= +c Ecriture NetCDF de l''entete du fichier "start_archive" +c======================================================================= + +c +c Preciser quelques attributs globaux: +c + ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 21, + & "Fichier start_archive") +c +c Definir les dimensions du fichiers: +c +c CHAMPS AJOUTES POUR LA VISUALISATION T,ps, etc... avec Grads ou ferret: + ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu) + ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv) + ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm) + ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx) + + ierr = NF_DEF_DIM (nid,"index", length, idim_index) + ierr = NF_DEF_DIM (nid,"rlonu", iip1, idim_rlonu) + ierr = NF_DEF_DIM (nid,"rlatv", jjm, idim_rlatv) + ierr = NF_DEF_DIM (nid,"interlayer", llmp1, idim_llmp1) + ierr = NF_DEF_DIM (nid,"Time", NF_UNLIMITED, idim_tim) + +c + ierr = NF_ENDDEF(nid) ! sortir du mode de definition + +c----------------------------------------------------------------------- +c Ecriture du tableau des parametres du run +c----------------------------------------------------------------------- + + call def_var(nid,"Time","Time","days since 00:00:00",1, + . idim_tim,nvarid,ierr) + + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"controle",NF_DOUBLE,1,idim_index,nvarid) +#else + ierr = NF_DEF_VAR (nid,"controle",NF_FLOAT,1,idim_index,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, + . "Parametres de controle") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl) +#endif + +c----------------------------------------------------------------------- +c Ecriture des longitudes et latitudes +c----------------------------------------------------------------------- + + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"rlonu",NF_DOUBLE,1,idim_rlonu,nvarid) +#else + ierr = NF_DEF_VAR (nid,"rlonu",NF_FLOAT,1,idim_rlonu,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, + . "Longitudes des points U") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu) +#endif +c + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"rlatu",NF_DOUBLE,1,idim_rlatu,nvarid) +#else + ierr = NF_DEF_VAR (nid,"rlatu",NF_FLOAT,1,idim_rlatu,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, + . "Latitudes des points U") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu) +#endif +c + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"rlonv",NF_DOUBLE,1,idim_rlonv,nvarid) +#else + ierr = NF_DEF_VAR (nid,"rlonv",NF_FLOAT,1,idim_rlonv,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, + . "Longitudes des points V") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv) +#endif +c + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"rlatv",NF_DOUBLE,1,idim_rlatv,nvarid) +#else + ierr = NF_DEF_VAR (nid,"rlatv",NF_FLOAT,1,idim_rlatv,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, + . "Latitudes des points V") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv) +#endif + +c----------------------------------------------------------------------- +c Ecriture des niveaux verticaux +c----------------------------------------------------------------------- + +c + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"ap",NF_DOUBLE,1,idim_llmp1,nvarid) +#else + ierr = NF_DEF_VAR (nid,"ap",NF_FLOAT,1,idim_llmp1,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 32, + . "Coef A: niveaux pression hybride") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,ap) +#endif +c + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"bp",NF_DOUBLE,1,idim_llmp1,nvarid) +#else + ierr = NF_DEF_VAR (nid,"bp",NF_FLOAT,1,idim_llmp1,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 35, + . "Coefficient B niveaux sigma hybride") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,bp) +#endif +c +c ---------------------- + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"aps",NF_DOUBLE,1,idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid,"aps",NF_FLOAT,1,idim_llm,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 36, + . "Coef AS: hybrid pressure in midlayers") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,aps) +#endif +c +c ---------------------- + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"bps",NF_DOUBLE,1,idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid,"bps",NF_FLOAT,1,idim_llm,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 30, + . "Coef BS: hybrid sigma midlayers") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,bps) +#endif +c +c ---------------------- + + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"presnivs",NF_DOUBLE,1,idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid,"presnivs",NF_FLOAT,1,idim_llm,nvarid) +#endif + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,presnivs) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,presnivs) +#endif +c ------------------------------------------------------------------ +c Variable uniquement pour visualisation avec Grads ou Ferret +c ------------------------------------------------------------------ + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"latitude",NF_DOUBLE,1,idim_rlatu,nvarid) +#else + ierr = NF_DEF_VAR (nid,"latitude",NF_FLOAT,1,idim_rlatu,nvarid) +#endif + ierr =NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14, + . "North latitude") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180) +#endif +c---------------------- + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr =NF_DEF_VAR(nid,"longitude", NF_DOUBLE, 1, idim_rlonv,nvarid) +#else + ierr = NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, idim_rlonv,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 14, + . "East longitude") + ierr = NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180) +#endif +c-------------------------- + ierr = NF_REDEF (nid) +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, + . idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, + . idim_llm,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",10,"pseudo-alt") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up") + + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt) +#endif + +!------------------------------- +! (soil) depth variable mlayer() (known from comsoil.h) +!------------------------------- + ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode + ! define variable +#ifdef NC_DOUBLE + ierr=NF_DEF_VAR(nid,"soildepth",NF_DOUBLE,1,idim_nsoilmx,nvarid) +#else + ierr=NF_DEF_VAR(nid,"soildepth",NF_FLOAT,1,idim_nsoilmx,nvarid) +#endif + ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20, + . "Soil mid-layer depth") + ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",1,"m") + ierr=NF_PUT_ATT_TEXT (nid,nvarid,"positive",4,"down") + ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode + ! write variable +#ifdef NC_DOUBLE + ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer) +#else + ierr=NF_PUT_VAR_REAL (nid,nvarid,mlayer) +#endif + +!--------------------- +! soil thermal inertia +!--------------------- + ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode + dims3(1)=idim_rlonv + dims3(2)=idim_rlatu + dims3(3)=idim_nsoilmx + ! define variable +#ifdef NC_DOUBLE + ierr=NF_DEF_VAR(nid,"inertiedat",NF_DOUBLE,3,dims3,nvarid) +#else + ierr=NF_DEF_VAR(nid,"inertiedat",NF_FLOAT,3,dims3,nvarid) +#endif + ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 20, + & "Soil thermal inertia") + ierr=NF_PUT_ATT_TEXT (nid,nvarid,"units",15, + & "J.s-1/2.m-2.K-1") + ierr=NF_ENDDEF(nid) ! Leave NetCDF define mode + ! write variable +#ifdef NC_DOUBLE + ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,ith) +#else + ierr=NF_PUT_VAR_REAL (nid,nvarid,ith) +#endif + +c----------------------------------------------------------------------- +c Ecriture aire et coefficients de passage cov. <-> contra. <--> naturel +c----------------------------------------------------------------------- + + ierr = NF_REDEF (nid) + dims2(1) = idim_rlonu + dims2(2) = idim_rlatu +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"cu",NF_DOUBLE,2,dims2,nvarid) +#else + ierr = NF_DEF_VAR (nid,"cu",NF_FLOAT,2,dims2,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, + . "Coefficient de passage pour U") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,cu) +#endif +c + ierr = NF_REDEF (nid) + dims2(1) = idim_rlonv + dims2(2) = idim_rlatv +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"cv",NF_DOUBLE,2,dims2,nvarid) +#else + ierr = NF_DEF_VAR (nid,"cv",NF_FLOAT,2,dims2,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, + . "Coefficient de passage pour V") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,cv) +#endif +c +c Aire de chaque maille: +c + ierr = NF_REDEF (nid) + dims2(1) = idim_rlonv + dims2(2) = idim_rlatu +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"aire",NF_DOUBLE,2,dims2,nvarid) +#else + ierr = NF_DEF_VAR (nid,"aire",NF_FLOAT,2,dims2,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, + . "Aires de chaque maille") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,aire) +#endif + +c----------------------------------------------------------------------- +c Ecriture du geopentiel au sol +c----------------------------------------------------------------------- + + ierr = NF_REDEF (nid) + dims2(1) = idim_rlonv + dims2(2) = idim_rlatu +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"phisinit",NF_DOUBLE,2,dims2,nvarid) +#else + ierr = NF_DEF_VAR (nid,"phisinit",NF_FLOAT,2,dims2,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19, + . "Geopotentiel au sol") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,phis) +#endif + + PRINT*,'iim,jjm,llm,idayref',iim,jjm,llm,idayref + PRINT*,'rad,omeg,g,mugaz,kappa', + s rad,omeg,g,43.49,kappa !mars temporaire (ecrire mugaz ensuite) + PRINT*,'daysec,dtvr',daysec,dtvr + + RETURN + END Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/iniphysiq_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/iniphysiq_mod.F90 (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/iniphysiq_mod.F90 (revision 1540) @@ -0,0 +1,189 @@ +MODULE iniphysiq_mod + +CONTAINS + +subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep, & + rlatu,rlatv,rlonu,rlonv,aire,cu,cv, & + prad,pg,pr,pcpp,iflag_phys) + +use dimphy, only : klev ! number of atmospheric levels +use mod_grid_phy_lmdz, only : klon_glo ! number of atmospheric columns + ! (on full grid) +use mod_phys_lmdz_para, only : klon_omp, & ! number of columns (on local omp grid) + klon_omp_begin, & ! start index of local omp subgrid + klon_omp_end, & ! end index of local omp subgrid + klon_mpi_begin ! start indes of columns (on local mpi grid) + +use comgeomphy, only : initcomgeomphy, & + airephy, & ! physics grid area (m2) + cuphy, & ! cu coeff. (u_covariant = cu * u) + cvphy, & ! cv coeff. (v_covariant = cv * v) + rlond, & ! longitudes + rlatd ! latitudes +use infotrac, only : nqtot ! number of advected tracers +use comgeomfi_h, only: ini_fillgeom +use temps_mod, only: day_ini, hour_ini +use phys_state_var_init_mod, only: phys_state_var_init +use regular_lonlat_mod, only: init_regular_lonlat, & + east, west, north, south, & + north_east, north_west, & + south_west, south_east + +implicit none + +include "iniprint.h" + +real,intent(in) :: prad ! radius of the planet (m) +real,intent(in) :: pg ! gravitational acceleration (m/s2) +real,intent(in) :: pr ! ! reduced gas constant R/mu +real,intent(in) :: pcpp ! specific heat Cp +real,intent(in) :: punjours ! length (in s) of a standard day +!integer,intent(in) :: ngrid ! number of horizontal grid points in the physics (full grid) +integer,intent(in) :: nlayer ! number of atmospheric layers +integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes +integer,intent(in) :: jj ! number of atompsheric columns along latitudes +real,intent(in) :: rlatu(jj+1) ! latitudes of the physics grid +real,intent(in) :: rlatv(jj) ! latitude boundaries of the physics grid +real,intent(in) :: rlonv(ii+1) ! longitudes of the physics grid +real,intent(in) :: rlonu(ii+1) ! longitude boundaries of the physics grid +real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2) +real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u) +real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v) +integer,intent(in) :: pdayref ! reference day of for the simulation +real,intent(in) :: ptimestep !physics time step (s) +integer,intent(in) :: iflag_phys ! type of physics to be called + +integer :: ibegin,iend,offset +integer :: i,j +character(len=20) :: modname='iniphysiq' +character(len=80) :: abort_message +real :: total_area_phy, total_area_dyn +real :: pi + +! boundaries, on global grid +real,allocatable :: boundslon_reg(:,:) +real,allocatable :: boundslat_reg(:,:) + +! global array, on full physics grid: +real,allocatable :: latfi(:) +real,allocatable :: lonfi(:) +real,allocatable :: cufi(:) +real,allocatable :: cvfi(:) +real,allocatable :: airefi(:) + +pi=2.*asin(1.0) + +IF (nlayer.NE.klev) THEN + write(*,*) 'STOP in ',trim(modname) + write(*,*) 'Problem with dimensions :' + write(*,*) 'nlayer = ',nlayer + write(*,*) 'klev = ',klev + abort_message = '' + CALL abort_gcm (modname,abort_message,1) +ENDIF + +!IF (ngrid.NE.klon_glo) THEN +! write(*,*) 'STOP in ',trim(modname) +! write(*,*) 'Problem with dimensions :' +! write(*,*) 'ngrid = ',ngrid +! write(*,*) 'klon = ',klon_glo +! abort_message = '' +! CALL abort_gcm (modname,abort_message,1) +!ENDIF + +! init regular global longitude-latitude grid points and boundaries +ALLOCATE(boundslon_reg(ii,2)) +ALLOCATE(boundslat_reg(jj+1,2)) + +DO i=1,ii + boundslon_reg(i,east)=rlonu(i) + boundslon_reg(i,west)=rlonu(i+1) +ENDDO + +boundslat_reg(1,north)= PI/2 +boundslat_reg(1,south)= rlatv(1) +DO j=2,jj + boundslat_reg(j,north)=rlatv(j-1) + boundslat_reg(j,south)=rlatv(j) +ENDDO +boundslat_reg(jj+1,north)= rlatv(jj) +boundslat_reg(jj+1,south)= -PI/2 + +! Write values in module regular_lonlat_mod +CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, & + boundslon_reg, boundslat_reg) + +! Generate global arrays on full physics grid +allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo)) +latfi(1)=rlatu(1) +lonfi(1)=0. +cufi(1) = cu(1) +cvfi(1) = cv(1) +DO j=2,jj + DO i=1,ii + latfi((j-2)*ii+1+i)= rlatu(j) + lonfi((j-2)*ii+1+i)= rlonv(i) + cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i) + cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i) + ENDDO +ENDDO +latfi(klon_glo)= rlatu(jj+1) +lonfi(klon_glo)= 0. +cufi(klon_glo) = cu((ii+1)*jj+1) +cvfi(klon_glo) = cv((ii+1)*jj-ii) + +! build airefi(), mesh area on physics grid +allocate(airefi(klon_glo)) +CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi) +! Poles are single points on physics grid +airefi(1)=sum(aire(1:ii,1)) +airefi(klon_glo)=sum(aire(1:ii,jj+1)) + +! Sanity check: do total planet area match between physics and dynamics? +total_area_dyn=sum(aire(1:ii,1:jj+1)) +total_area_phy=sum(airefi(1:klon_glo)) +IF (total_area_dyn/=total_area_phy) THEN + WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!' + WRITE (lunout, *) ' in the dynamics total_area_dyn=', total_area_dyn + WRITE (lunout, *) ' but in the physics total_area_phy=', total_area_phy + IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN + ! stop here if the relative difference is more than 0.001% + abort_message = 'planet total surface discrepancy' + CALL abort_gcm(modname, abort_message, 1) + ENDIF +ENDIF + + + +!$OMP PARALLEL +! Now generate local lon/lat/cu/cv/area arrays +call initcomgeomphy + +!!!!$OMP PARALLEL PRIVATE(ibegin,iend) +!!!$OMP+ SHARED(parea,pcu,pcv,plon,plat) + +offset=klon_mpi_begin-1 +airephy(1:klon_omp)=airefi(offset+klon_omp_begin:offset+klon_omp_end) +cuphy(1:klon_omp)=cufi(offset+klon_omp_begin:offset+klon_omp_end) +cvphy(1:klon_omp)=cvfi(offset+klon_omp_begin:offset+klon_omp_end) +rlond(1:klon_omp)=lonfi(offset+klon_omp_begin:offset+klon_omp_end) +rlatd(1:klon_omp)=latfi(offset+klon_omp_begin:offset+klon_omp_end) + +! copy some fundamental parameters to physics +! and do some initializations +call phys_state_var_init(klon_omp,nlayer,nqtot, & + day_ini,hour_ini,punjours,ptimestep, & + prad,pg,pr,pcpp) +call ini_fillgeom(klon_omp,rlatd,rlond,airephy) +call conf_phys(klon_omp,nlayer,nqtot) + +! Initialize some "temporal and calendar" related variables +!CALL init_time(day_ini,hour_ini,punjours,ptimestep) + +!$OMP END PARALLEL + + +end subroutine iniphysiq + + +END MODULE iniphysiq_mod Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/interp_vert.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/interp_vert.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/interp_vert.F (revision 1540) @@ -0,0 +1,70 @@ +c****************************************************** + SUBROUTINE interp_vert(varo,varn,lmo,lmn,apso,bpso, + & aps,bps,ps,Nhoriz) +c +c interpolation lineaire pour passer +c a une nouvelle discretisation verticale pour +c les variables de GCM +c Francois Forget (01/1995) +c Modif pour coordonnees hybrides FF (03/2003) +c********************************************************** + + IMPLICIT NONE + +c Declarations: +c ============== +c +c ARGUMENTS +c """"""""" + + integer lmo ! dimensions ancienne couches (input) + integer lmn ! dimensions nouvelle couches (input) + + real apso(lmo),bpso(lmo)! anciennes coord hybrides midlayer (input) + real aps(lmn), bps(lmn)! nouvelles coord hybrides (midlayer) (input) + + integer Nhoriz ! nombre de point horizontale (input) + real ps(nhoriz) !pression de surface (input) + + real varo(Nhoriz,lmo) ! var dans l''ancienne grille (input) + real varn(Nhoriz,lmn) ! var dans la nouvelle grille (output) + +c Autres variables +c """""""""""""""" + integer n, ln ,lo + real coef + REAL sigmo(lmo) ! niveau sigma des variables dans les anciennes coord + REAL sigmn(lmn) ! niveau sigma des variables dans les nouvelles coord + +c run +c ==== + + do n=1,Nhoriz + do ln=1,lmn + sigmn(ln)=aps(ln)/ps(n)+bps(ln) + end do + do lo=1,lmo + sigmo(lo)=apso(lo)/ps(n)+bpso(lo) + end do + + do ln=1,lmn + if (sigmn(ln).ge.sigmo(1))then + varn(n,ln) = varo(n,1) + else if (sigmn(ln).le.sigmo(lmo)) then + varn(n,ln) = varo(n,lmo) + else + do lo =1,lmo-1 + if ( (sigmn(ln).le.sigmo(lo)).and. + & (sigmn(ln).gt.sigmo(lo+1)) )then + coef = (sigmn(ln)-sigmo(lo))/(sigmo(lo+1)-sigmo(lo)) + varn(n,ln)=varo(n,lo) +coef*(varo(n,lo+1)-varo(n,lo)) + end if + end do + end if + end do + + end do + + + return + end Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F =================================================================== --- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F (revision 1540) +++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F (revision 1540) @@ -0,0 +1,1428 @@ + SUBROUTINE lect_start_archive(ngrid,nlayer,nqtot, + & date,tsurf,tsoil,emis,q2, + & t,ucov,vcov,ps,co2ice,h,phisold_newgrid, + & q,qsurf,tauscaling,surfith,nid) +c======================================================================= +c +c +c Auteur: 05/1997 , 12/2003 : coord hybride FF +c ------ +c +c +c Objet: Lecture des variables d'un fichier "start_archive" +c Plus besoin de régler ancienne valeurs grace +c a l'allocation dynamique de memoire (Yann Wanherdrick) +c +c +c +c======================================================================= + use infotrac, only: tname + use comsoil_h, only: nsoilmx, layer, mlayer, volcapa, inertiedat + use planete_h + USE comvert_mod, ONLY: ap,bp,aps,bps,preff + USE comconst_mod, ONLY: kappa,g,pi + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom2.h" +#include "netcdf.inc" +c======================================================================= +c Declarations +c======================================================================= + +! routine arguments +! ----------------- + integer,intent(in) :: ngrid ! number of atmosphferic columns + ! on new physics grid + integer,intent(in) :: nlayer ! number of atmospheric layers + ! on new grid + integer,intent(in) :: nqtot ! number of advected tracers + REAL,INTENT(OUT) :: date + REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants + REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1) + REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) + REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature + REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature + REAL,INTENT(OUT) :: co2ice(ngrid) ! CO2 ice layer + REAL,INTENT(OUT) :: emis(ngrid) + REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot) + REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor + REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1) + REAL,INTENT(OUT) :: t(iip1,jjp1,llm) + real,intent(in) :: surfith(iip1,jjp1) ! surface thermal inertia + INTEGER,INTENT(IN) :: nid ! input NetCDF file identifier + + + +c Old variables dimensions (from file) +c------------------------------------ + INTEGER imold,jmold,lmold,nsoilold,nqold + +c Variables pour les lectures des fichiers "ini" +c-------------------------------------------------- +! INTEGER sizei, + integer timelen,dimid +! INTEGER length +! parameter (length = 100) + INTEGER tab0 + INTEGER isoil,iq,iqmax + CHARACTER*2 str2 + +! REAL dimfirst(4) ! tableau contenant les 1ers elements des dimensions + +! REAL dimlast(4) ! tableau contenant les derniers elements des dimensions + +! REAL dimcycl(4) ! tableau contenant les periodes des dimensions +! CHARACTER*120 dimsource +! CHARACTER*16 dimname +! CHARACTER*80 dimtitle +! CHARACTER*40 dimunits +! INTEGER dimtype + +! INTEGER dimord(4) ! tableau contenant l''ordre +! data dimord /1,2,3,4/ ! de sortie des dimensions + +! INTEGER vardim(4) + INTEGER memo +! character (len=50) :: tmpname + +c Variable histoire +c------------------ + REAL ::qtot(iip1,jjp1,llm) + +c autre variables dynamique nouvelle grille +c------------------------------------------ + +c!-*- +! integer klatdat,klongdat +! PARAMETER (klatdat=180,klongdat=360) + +c Physique sur grille scalaire +c---------------------------- + +c variable physique +c------------------ +c REAL phisfi(ngrid) + + INTEGER i,j,l + INTEGER nvarid +c REAL year_day,periheli,aphelie,peri_day +c REAL obliquit,z0,emin_turb,lmixmin +c REAL emissiv,emisice(2),albedice(2),tauvis +c REAL iceradius(2) , dtemisice(2) + +! EXTERNAL RAN1 +! REAL RAN1 +! EXTERNAL geopot,inigeom + integer ierr +! integer ismin +! external ismin +! CHARACTER*80 datapath + integer, dimension(4) :: start,count + +c Variable nouvelle grille naturelle au point scalaire +c------------------------------------------------------ + real us(iip1,jjp1,llm),vs(iip1,jjp1,llm) + real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx) + real inertiedatS(iip1,jjp1,nsoilmx) + real co2iceS(iip1,jjp1),emisS(iip1,jjp1) + REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot) + real tauscalingS(iip1,jjp1) + + real ptotal, co2icetotal + +c Var intermediaires : vent naturel, mais pas coord scalaire +c----------------------------------------------------------- + real vnat(iip1,jjm,llm),unat(iip1,jjp1,llm) + + +c Variable de l'ancienne grille +c--------------------------------------------------------- + + real, dimension(:), allocatable :: timelist + real, dimension(:), allocatable :: rlonuold, rlatvold + real, dimension(:), allocatable :: rlonvold, rlatuold + real, dimension(:), allocatable :: apsold,bpsold + real, dimension(:), allocatable :: mlayerold + real, dimension(:,:,:), allocatable :: uold,vold,told,q2old + real, dimension(:,:,:), allocatable :: tsoilold,qsurfold + real, dimension(:,:,:),allocatable :: tsoiloldnew +! tsoiloldnew: old soil values, but along new subterranean grid + real, dimension(:,:,:), allocatable :: inertiedatold +! inertiedatoldnew: old inertia values, but along new subterranean grid + real, dimension(:,:,:), allocatable :: inertiedatoldnew + real, dimension(:,:), allocatable :: psold,phisold + real, dimension(:,:), allocatable :: co2iceold,tsurfold + real, dimension(:,:), allocatable :: emisold + real, dimension(:,:,:,:), allocatable :: qold + real, dimension(:,:), allocatable :: tauscalingold + + real tab_cntrl(100) + + real ptotalold, co2icetotalold + + logical :: olddepthdef=.false. ! flag +! olddepthdef=.true. if soil depths are in 'old' (unspecified) format + logical :: depthinterpol=.false. ! flag +! depthinterpol=.true. if interpolation will be requiered + logical :: therminertia_3D=.true. ! flag +! therminertia_3D=.true. if thermal inertia is 3D and read from datafile +c Variable intermediaires iutilise pour l'extrapolation verticale +c---------------------------------------------------------------- + real, dimension(:,:,:), allocatable :: var,varp1 + real, dimension(:), allocatable :: oldgrid, oldval + real, dimension(:), allocatable :: newval +! real, dimension(:), allocatable :: oldmlayer + +! real surfithfi(ngrid) + ! surface thermal inertia at old horizontal grid resolution + real, dimension(:,:), allocatable :: surfithold + +! flag which identifies if archive file is using old tracer names (qsurf01,...) + logical :: oldtracernames=.false. + integer :: counter + character(len=30) :: txt ! to store some text + real :: tmpval ! to store a temporary variable/value + +c======================================================================= + +! 0. Preliminary stuff + +! check if tracers follow old naming convention (q01, q02, q03, ...) + counter=0 + do iq=1,nqtot + txt= " " + write(txt,'(a1,i2.2)')'q',iq + ierr=NF_INQ_VARID(nid,txt,nvarid) + if (ierr.ne.NF_NOERR) then + ! did not find old tracer name + exit ! might as well stop here + else + ! found old tracer name + counter=counter+1 + endif + enddo + if (counter.eq.nqtot) then + write(*,*) "lect_start_archive: tracers seem to follow old ", + & "naming convention (q01, q02,...)" + oldtracernames=.true. + endif + + +!----------------------------------------------------------------------- +! 1. Read data dimensions (i.e. size and length) +!----------------------------------------------------------------------- + +! 1.2 Read the various dimension lengths of data in file + + ierr= NF_INQ_DIMID(nid,"Time",dimid) + if (ierr.ne.NF_NOERR) then + ierr= NF_INQ_DIMID(nid,"temps",dimid) + endif + ierr= NF_INQ_DIMLEN(nid,dimid,timelen) + if (ierr.ne.NF_NOERR) then + write(*,*) 'lect_start_archive error: cannot find Time length' + stop + else + write(*,*) "lect_start_archive: timelen=",timelen + endif + + ierr= NF_INQ_DIMID(nid,"latitude",dimid) + if (ierr.ne.NF_NOERR) then + ierr= NF_INQ_DIMID(nid,"rlatu",dimid) + endif + ierr= NF_INQ_DIMLEN(nid,dimid,jmold) + if (ierr.ne.NF_NOERR) then + write(*,*) 'lect_start_archive error: cannot find lat length' + stop + else + write(*,*) "lect_start_archive: jmold=",jmold + endif + jmold=jmold-1 + + ierr= NF_INQ_DIMID(nid,"longitude",dimid) + if (ierr.ne.NF_NOERR) then + ierr= NF_INQ_DIMID(nid,"rlonv",dimid) + endif + ierr= NF_INQ_DIMLEN(nid,dimid,imold) + if (ierr.ne.NF_NOERR) then + write(*,*) 'lect_start_archive error: cannot find lon length' + stop + else + write(*,*) "lect_start_archive: imold=",imold + endif + imold=imold-1 + + ierr= NF_INQ_DIMID(nid,"altitude",dimid) + if (ierr.ne.NF_NOERR) then + ierr= NF_INQ_DIMID(nid,"sig_s",dimid) + endif + ierr= NF_INQ_DIMLEN(nid,dimid,lmold) + if (ierr.ne.NF_NOERR) then + write(*,*) 'lect_start_archive error: cannot find alt length' + stop + else + write(*,*) "lect_start_archive: lmold=",lmold + endif + + nqold=0 + do + write(str2,'(i2.2)') nqold+1 + ierr= NF_INQ_VARID(nid,'q'//str2,dimid) +! write(*,*) 'q'//str2 + if (ierr.eq.NF_NOERR) then + nqold=nqold+1 + else + exit + endif + enddo + +! 1.2.1 find out the # of subsurface_layers + nsoilold=0 !dummy initialisation + ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid) + if (ierr.eq.NF_NOERR) then + ierr=NF_INQ_DIMLEN(nid,dimid,nsoilold) + if (ierr.ne.NF_NOERR) then + write(*,*)'lec_start_archive: ', + & 'Failed reading subsurface_layers length' + endif + else + write(*,*)"lec_start_archive: did not find subsurface_layers" + endif + + if (nsoilold.eq.0) then ! 'old' archive format; + ! must use Tg//str2 fields to compute nsoilold + write(*,*)"lec_start_archive: building nsoilold from Tg fields" + do + write(str2,'(i2.2)') nsoilold+1 + ierr=NF_INQ_VARID(nid,'Tg'//str2,dimid) + if (ierr.eq.NF_NOERR) then + nsoilold=nsoilold+1 + else + exit + endif + enddo + endif + + + if (nsoilold.ne.nsoilmx) then ! interpolation will be required + depthinterpol=.true. + endif + +! 1.3 Report dimensions + + write(*,*) "lect_start_archive: Start_archive dimensions:" + write(*,*) "longitude: ",imold + write(*,*) "latitude: ",jmold + write(*,*) "altitude: ",lmold + if (nqold.gt.0) then + write(*,*) "old tracers q*: ",nqold + endif + write(*,*) "subsurface_layers: ",nsoilold + if (depthinterpol) then + write(*,*) " => Warning, nsoilmx= ",nsoilmx + write(*,*) ' which implies that you want subterranean interpola + &tion.' + write(*,*) ' Otherwise, set nsoilmx -in comsoil_h- to: ',nsoilold + endif + write(*,*) "time lenght: ",timelen + write(*,*) + +!----------------------------------------------------------------------- +! 2. Allocate arrays to store datasets +!----------------------------------------------------------------------- + + allocate(timelist(timelen)) + allocate(rlonuold(imold+1), rlatvold(jmold)) + allocate(rlonvold(imold+1), rlatuold(jmold+1)) + allocate (apsold(lmold),bpsold(lmold)) + allocate(uold(imold+1,jmold+1,lmold)) + allocate(vold(imold+1,jmold+1,lmold)) + allocate(told(imold+1,jmold+1,lmold)) + allocate(psold(imold+1,jmold+1)) + allocate(phisold(imold+1,jmold+1)) + allocate(qold(imold+1,jmold+1,lmold,nqtot)) + allocate(co2iceold(imold+1,jmold+1)) + allocate(tsurfold(imold+1,jmold+1)) + allocate(emisold(imold+1,jmold+1)) + allocate(q2old(imold+1,jmold+1,lmold+1)) +! allocate(tsoilold(imold+1,jmold+1,nsoilmx)) + allocate(tsoilold(imold+1,jmold+1,nsoilold)) + allocate(tsoiloldnew(imold+1,jmold+1,nsoilmx)) + allocate(inertiedatold(imold+1,jmold+1,nsoilold)) ! soil thermal inertia + allocate(inertiedatoldnew(imold+1,jmold+1,nsoilmx)) + ! surface thermal inertia at old horizontal grid resolution + allocate(surfithold(imold+1,jmold+1)) + allocate(mlayerold(nsoilold)) + allocate(qsurfold(imold+1,jmold+1,nqtot)) + allocate(tauscalingold(imold+1,jmold+1)) + + allocate(var (imold+1,jmold+1,llm)) + allocate(varp1 (imold+1,jmold+1,llm+1)) + + write(*,*) 'q2',ngrid,nlayer+1 + write(*,*) 'q2S',iip1,jjp1,llm+1 + write(*,*) 'q2old',imold+1,jmold+1,lmold+1 + +!----------------------------------------------------------------------- +! 3. Read time-independent data +!----------------------------------------------------------------------- + +C----------------------------------------------------------------------- +c 3.1. Lecture du tableau des parametres du run +c (pour la lecture ulterieure de "ptotalold" et "co2icetotalold") +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "controle", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "Lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF +c + tab0 = 50 + +c----------------------------------------------------------------------- +c 3.2 Lecture des longitudes et latitudes +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "rlonv", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonvold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlonvold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlatu", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatuold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlatuold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlonu", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonuold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlonuold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF +c + ierr = NF_INQ_VARID (nid, "rlatv", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatvold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, rlatvold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF +c + +c----------------------------------------------------------------------- +c 3.3. Lecture des niveaux verticaux +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "aps", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + apsold=0 + PRINT*, " set to 0" + ELSE +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, apsold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, apsold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + ENDIF + ENDIF +c + ierr = NF_INQ_VARID (nid, "bps", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + PRINT*, "It must be an old start_archive, lets look for sig_s" + ierr = NF_INQ_VARID (nid, "sig_s", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "Nothing to do..." + CALL abort + ENDIF + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bpsold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, bpsold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + END IF + +c----------------------------------------------------------------------- +c 3.4 Read Soil layers depths +c----------------------------------------------------------------------- + + ierr=NF_INQ_VARID(nid,"soildepth",nvarid) + if (ierr.ne.NF_NOERR) then + write(*,*)'lect_start_archive: Could not find ' + write(*,*)' => Assuming this is an archive in old format' + olddepthdef=.true. + depthinterpol=.true. + ! this is how soil depth was defined in ye old days + do isoil=1,nsoilold + mlayerold(isoil)=sqrt(887.75/3.14)*((2.**(isoil-0.5))-1.) + enddo + else +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayerold) +#else + ierr = NF_GET_VAR_REAL(nid,nvarid,mlayerold) +#endif + if (ierr .NE. NF_NOERR) then + PRINT*, "lect_start_archive: Failed reading " + CALL abort + endif + + endif !of if(ierr.ne.NF_NOERR) + + ! Read (or build) mlayer() + if (depthinterpol) then + ! Build (default) new soil depths (mlayer(:) is in comsoil.h), + ! as in soil_settings.F + write(*,*)' => Building default soil depths' + do isoil=0,nsoilmx-1 + mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) + enddo + write(*,*)' => mlayer: ',mlayer + ! Also build (default) new soil interlayer depth layer(:) + do isoil=1,nsoilmx + layer(isoil)=sqrt(mlayer(0)*mlayer(1))* + & ((mlayer(1)/mlayer(0))**(isoil-1)) + enddo + write(*,*)' => layer: ',layer + else ! read mlayer() from file +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer) +#else + ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer) +#endif + if (ierr .NE. NF_NOERR) then + PRINT*, "lect_start_archive: Failed reading " + CALL abort + endif + ! Also build (default) soil interlayer depth layer(:) + do isoil=1,nsoilmx + layer(isoil)=sqrt(mlayer(0)*mlayer(1))* + & ((mlayer(1)/mlayer(0))**(isoil-1)) + enddo + endif ! of if (depthinterpol) + +c----------------------------------------------------------------------- +c 3.5 Read Soil thermal inertia +c----------------------------------------------------------------------- + + ierr=NF_INQ_VARID(nid,"inertiedat",nvarid) + if (ierr.ne.NF_NOERR) then + write(*,*)'lect_start_archive: Could not find ' + write(*,*)' => Assuming this is an archive in old format' + therminertia_3D=.false. + write(*,*)' => Thermal inertia will be read from reference file' + volcapa=1.e6 + write(*,*)' and soil volumetric heat capacity is set to ', + & volcapa + else +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedatold) +#else + ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedatold) +#endif + if (ierr .NE. NF_NOERR) then + PRINT*, "lect_start_archive: Failed reading " + CALL abort + endif + endif + +c----------------------------------------------------------------------- +c 3.6 Lecture geopotentiel au sol +c----------------------------------------------------------------------- +c + ierr = NF_INQ_VARID (nid, "phisinit", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: is missing" + CALL abort + ENDIF +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisold) +#else + ierr = NF_GET_VAR_REAL(nid, nvarid, phisold) +#endif + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: Failed loading " + CALL abort + ENDIF + +C----------------------------------------------------------------------- +c lecture de "ptotalold" et "co2icetotalold" +c----------------------------------------------------------------------- + ptotalold = tab_cntrl(tab0+49) + co2icetotalold = tab_cntrl(tab0+50) + +c----------------------------------------------------------------------- +c 4. Lecture du temps et choix +c----------------------------------------------------------------------- + +c lecture du temps +c + ierr = NF_INQ_DIMID (nid, "Time", nvarid) + IF (ierr .NE. NF_NOERR) THEN + ierr = NF_INQ_DIMID (nid, "temps", nvarid) + IF (ierr .NE. NF_NOERR) THEN + PRINT*, "lect_start_archive: