Index: trunk/LMDZ.MARS/libf/phymars/testphys1d.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 890) +++ trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 899) @@ -29,4 +29,5 @@ #include "comgeomfi.h" #include "surfdat.h" +#include "slope.h" #include "comsoil.h" #include "comdiurn.h" @@ -100,46 +101,46 @@ c ------------------------------------------------------ -c Constantes prescrites ICI +c Prescribed constants to be set here c ------------------------------------------------------ pi=2.E+0*asin(1.E+0) -c Constante de la Planete Mars +c Mars planetary constants c ---------------------------- - rad=3397200. ! rayon de mars (m) ~3397200 m - daysec=88775. ! duree du sol (s) ~88775 s - omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1) - g=3.72 ! gravite (m.s-2) ~3.72 - mugaz=43.49 ! Masse molaire de l'atm (g.mol-1) ~43.49 + rad=3397200. ! mars radius (m) ~3397200 m + daysec=88775. ! length of a sol (s) ~88775 s + omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1) + g=3.72 ! gravity (m.s-2) ~3.72 + mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49 rcp=.256793 ! = r/cp ~0.256793 r= 8.314511E+0 *1000.E+0/mugaz cpp= r/rcp - year_day = 669 ! duree de l'annee (sols) ~668.6 - periheli = 206.66 ! dist.min. soleil-mars (Mkm) ~206.66 - aphelie = 249.22 ! dist.max. soleil-mars (Mkm) ~249.22 - peri_day = 485. ! date du perihelie (sols depuis printemps) - obliquit = 25.2 ! Obliquite de la planete (deg) ~25.2 - -c Parametres Couche limite et Turbulence -c -------------------------------------- - z0_default = 1.e-2 ! surface roughness (m) ~0.01 - emin_turb = 1.e-6 ! energie minimale ~1.e-8 - lmixmin = 30 ! longueur de melange ~100 - -c propriete optiques des calottes et emissivite du sol + year_day = 669 ! lenght of year (sols) ~668.6 + periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66 + aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22 + peri_day = 485. ! perihelion date (sols since N. Spring) + obliquit = 25.2 ! Obliquity (deg) ~25.2 + +c Planetary Boundary Layer and Turbulence parameters +c -------------------------------------------------- + z0_default = 1.e-2 ! surface roughness (m) ~0.01 + emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8 + lmixmin = 30 ! mixing length ~100 + +c cap properties and surface emissivities c ---------------------------------------------------- - emissiv= 0.95 ! Emissivite du sol martien ~.95 - emisice(1)=0.95 ! Emissivite calotte nord - emisice(2)=0.95 ! Emissivite calotte sud - albedice(1)=0.5 ! Albedo calotte nord - albedice(2)=0.5 ! Albedo calotte sud + emissiv= 0.95 ! Bare ground emissivity ~.95 + emisice(1)=0.95 ! Northern cap emissivity + emisice(2)=0.95 ! Southern cap emisssivity + albedice(1)=0.5 ! Northern cap albedo + albedice(2)=0.5 ! Southern cap albedo iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north) iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south) dtemisice(1) = 2. ! time scale for snow metamorphism (north) dtemisice(2) = 2. ! time scale for snow metamorphism (south - hybrid=.false. + c ------------------------------------------------------ -c Lecture des parametres dans "run.def" +c Loading run parameters from "run.def" file c ------------------------------------------------------ @@ -354,4 +355,7 @@ tnom(iq)=str7 enddo + ! and just to be clean, also initialize tracers to zero for physdem1 + q(:,:)=0 + qsurf(:)=0 endif ! of if (tracer) @@ -359,7 +363,7 @@ !write(*,*) "testphys1d qsurf", qsurf -c Date et heure locale du debut du run -c ------------------------------------ -c Date (en sols depuis le solstice de printemps) du debut du run +c Date and local time at beginning of run +c --------------------------------------- +c Date (in sols since spring solstice) at beginning of run day0 = 0 ! default value for day0 write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?' @@ -367,5 +371,5 @@ day=float(day0) write(*,*) " day0 = ",day0 -c Heure de demarrage +c Local time at beginning of run time=0 ! default value for time write(*,*)'Initial local time (in hours, between 0 and 24)?' @@ -374,5 +378,5 @@ time=time/24.E+0 ! convert time (hours) to fraction of sol -c Discretisation (Definition de la grille et des pas de temps) +c Discretization (Definition of grid and time steps) c -------------- c @@ -393,5 +397,6 @@ ndt=ndt*day_step dtphys=daysec/day_step -c Pression de surface sur la planete + +c Imposed surface pressure c ------------------------------------ c @@ -400,11 +405,11 @@ call getin("psurf",psurf) write(*,*) " psurf = ",psurf -c Pression de reference - pa=20. - preff=610. - -c Proprietes des poussiere aerosol +c Reference pressures + pa=20. ! transition pressure (for hybrid coord.) + preff=610. ! reference surface pressure + +c Aerosol properties c -------------------------------- - tauvis=0.2 ! default value for tauvis + tauvis=0.2 ! default value for tauvis (dust opacity) write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref call getin("tauvis",tauvis) @@ -439,6 +444,6 @@ long(1)=long(1)*pi/180.E+0 -c Initialisation albedo / inertie du sol -c -------------------------------------- +c Initialize albedo / soil thermal inertia +c ---------------------------------------- c albedodat(1)=0.2 ! default value for albedodat @@ -456,6 +461,17 @@ call getin("z0",z0(1)) write(*,*) " z0 = ",z0(1) -c -c pour le schema d'ondes de gravite + +! Initialize local slope parameters (only matters if "callslope" +! is .true. in callphys.def) + ! slope inclination angle (deg) 0: horizontal, 90: vertical + theta_sl(1)=0.0 ! default: no inclination + call getin("slope_inclination",theta_sl(1)) + ! slope orientation (deg) + ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward + psi_sl(1)=0.0 ! default value + call getin("slope_orientation",psi_sl(1)) + +c +c for the gravity wave scheme c --------------------------------- c @@ -467,27 +483,25 @@ - -c Initialisation speciales "physiq" -c --------------------------------- -c la surface de chaque maille est inutile en 1D ---> +c Specific initializations for "physiq" +c ------------------------------------- +c mesh surface (not a very usefull quantity in 1D) area(1)=1.E+0 -c le geopotentiel au sol est inutile en 1D car tout est controle -c par la pression de surface ---> +c surface geopotential is not used (or useful) since in 1D +c everything is controled by surface pressure phisfi(1)=0.E+0 -c "inifis" reproduit un certain nombre d'initialisations deja faites -c + lecture des clefs de callphys.def -c + calcul de la frequence d'appel au rayonnement -c + calcul des sinus et cosinus des longitude latitude +c "inifis" does some initializations (some of which have already been +c done above!) and loads parameters set in callphys.def !Mars possible matter with dtphys in input and include!!! CALL inifis(1,llm,day0,daysec,dtphys, . lati,long,area,rad,g,r,cpp) -c Initialisation pour prendre en compte les vents en 1-D + +c Initialization to take into account prescribed winds c ------------------------------------------------------ ptif=2.E+0*omeg*sinlat(1) -c vent geostrophique +c geostrophic wind gru=10. ! default value for gru PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?' @@ -500,5 +514,5 @@ write(*,*) " v = ",grv -c Initialisation des vents au premier pas de temps +c Initialize winds for first time step DO ilayer=1,nlayer u(ilayer)=gru @@ -506,10 +520,10 @@ ENDDO -c energie cinetique turbulente +c Initialize turbulente kinetic energy DO ilevel=1,nlevel q2(ilevel)=0.E+0 ENDDO -c glace de CO2 au sol +c CO2 ice on the surface c ------------------- co2ice=0.E+0 ! default value for co2ice @@ -519,5 +533,5 @@ c -c emissivite +c emissivity c ---------- emis=emissiv @@ -529,5 +543,5 @@ -c calcul des pressions et altitudes en utilisant les niveaux sigma +c Compute pressures and altitudes of atmospheric levels c ---------------------------------------------------------------- @@ -555,9 +569,9 @@ -c profil de temperature au premier appel +c Initialize temperature profile c -------------------------------------- pks=psurf**rcp -c altitude en km dans profile: on divise zlay par 1000 +c altitude in km in profile: divide zlay by 1000 tmp1(0)=0.E+0 DO ilayer=1,nlayer @@ -591,10 +605,10 @@ enddo -c Initialisation des traceurs +c Initialize traceurs c --------------------------- if (photochem.or.callthermos) then - write(*,*) 'Especes chimiques initialisees' - ! thermo=0: initialisation pour toutes les couches + write(*,*) 'Initializing chemical species' + ! thermo=0: initialize over all atmospheric layers thermo=0 call inichim_newstart(q,psurf,sig,nqmx,lati,long,area, @@ -602,7 +616,7 @@ endif -c Regarder si le sol est un reservoir de glace d'eau +c Check if the surface is a water ice reservoir c -------------------------------------------------- - watercaptag(ngridmx)=.false. ! Par defaut il n'y pas de glace au sol + watercaptag(ngridmx)=.false. ! Default: no water ice reservoir print *,'Water ice cap on ground ?' call getin("watercaptag",watercaptag) @@ -610,8 +624,7 @@ -c Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl" +c Initialization for GRADS outputs in "g1d.dat" and "g1d.ctl" c ---------------------------------------------------------------- -c (effectuee avec "writeg1d" a partir de "physiq.F" -c ou tout autre subroutine physique, y compris celle ci). +c (output done in "writeg1d", typically called by "physiq.F") g1d_nlayer=nlayer @@ -623,9 +636,8 @@ g2d_premier=.true. -c Ecriture de "startfi" +c Write a "startfi" file c -------------------- -c (Ce fichier sera aussitot relu au premier -c appel de "physiq", mais il est necessaire pour passer -c les variables purement physiques a "physiq"... +c This file will be read during the first call to "physiq". +c It is needed to transfert physics variables to "physiq"... call physdem1("startfi.nc",long,lati,nsoilmx,nqmx, @@ -633,6 +645,7 @@ . tsoil,co2ice,emis,q2,qsurf,area,albedodat, . inertiedat,zmea,zstd,zsig,zgam,zthe) + c======================================================================= -c BOUCLE TEMPORELLE DU MODELE 1D +c 1D MODEL TIME STEPPING LOOP c======================================================================= c @@ -656,5 +669,5 @@ ENDIF -c calcul du geopotentiel +c compute geopotential c ~~~~~~~~~~~~~~~~~~~~~ DO ilayer=1,nlayer @@ -670,5 +683,5 @@ ENDDO -c appel de la physique +c call physics c -------------------- ! write(*,*) "testphys1d avant q", q(1,:) @@ -679,14 +692,14 @@ , u, v,temp, q, , w, -C - sorties +C - outputs s du, dv, dtemp, dq,dpsurf,tracerdyn) ! write(*,*) "testphys1d apres q", q(1,:) -c evolution du vent : modele 1D -c ----------------------------- - -c la physique calcule les derivees temporelles de u et v. -c on y rajoute betement un effet Coriolis. +c wind increment : specific for 1D +c -------------------------------- + +c The physics compute the tendencies on u and v, +c here we just add Coriolos effect c c DO ilayer=1,nlayer @@ -695,5 +708,5 @@ c ENDDO -c Pour certain test : pas de coriolis a l'equateur +c For some tests : No coriolis force at equator c if(lati(1).eq.0.) then DO ilayer=1,nlayer @@ -704,5 +717,5 @@ c c -c Calcul du temps au pas de temps suivant +c Compute time for next time step c --------------------------------------- firstcall=.false. @@ -713,5 +726,5 @@ ENDIF -c calcul des vitesses et temperature au pas de temps suivant +c compute winds and temperature for next time step c ---------------------------------------------------------- @@ -722,8 +735,8 @@ ENDDO -c calcul des pressions au pas de temps suivant +c compute pressure for next time step c ---------------------------------------------------------- - psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface + psurf=psurf+dtphys*dpsurf ! surface pressure change DO ilevel=1,nlevel plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) @@ -733,5 +746,5 @@ ENDDO -! increment tracer values +! increment tracers DO iq = 1, nqmx DO ilayer=1,nlayer @@ -740,13 +753,11 @@ ENDDO - ENDDO ! fin de la boucle temporelle + ENDDO ! of idt=1,ndt ! end of time stepping loop c ======================================================== -c GESTION DES SORTIE +c OUTPUTS c ======================================================== -c fermeture pour conclure les sorties format grads dans "g1d.dat" -c et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F" -c ou tout autre subroutine physique, y compris celle ci). +c finalize and close grads files "g1d.dat" and "g1d.ctl" c CALL endg1d(1,nlayer,zphi/(g*1000.),ndt) @@ -758,10 +769,23 @@ c*********************************************************************** c*********************************************************************** -c Subroutines Bidons utilise seulement en 3D, mais -c necessaire a la compilation de testphys1d en 1-D - - subroutine gr_fi_dyn - RETURN - END +c Dummy subroutines used only in 3D, but required to +c compile testphys1d (to cleanly use writediagfi) + + subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn) + + IMPLICIT NONE + + INTEGER im,jm,ngrid,nfield + REAL pdyn(im,jm,nfield) + REAL pfi(ngrid,nfield) + + if (ngrid.ne.1) then + write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!" + stop + endif + + pdyn(1,1,1:nfield)=pfi(1,1:nfield) + + end c***********************************************************************