Changeset 705 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Jun 14, 2012, 4:14:11 PM (13 years ago)

Author:

emillour

Message:

Mars GCM:

• Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
• Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
• Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• callkeys.h (modified) (2 diffs)
• inifis.F (modified) (2 diffs)
• physiq.F (modified) (7 diffs)

trunk/LMDZ.MARS/libf/phymars/callkeys.h

@@ -16 +16 @@
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
-     &   ,dustbin,nltemodel,nircorr
+     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear
      
       COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff,     &
@@ -47 +47 @@
       integer ilwn
       integer ncouche
+      integer solvarmod   ! model for solar EUV variation
+      integer solvaryear  ! mars year for realisticly varying solar EUV 
 
       logical rayleigh

trunk/LMDZ.MARS/libf/phymars/inifis.F

@@ -584 +584 @@
          write(*,*) " thermochem = ",thermochem
 
+         write(*,*) "Method to include solar variability"
+         write(*,*) "0-> old method (using solarcondate); ",
+     &                  "1-> variability wit E10.7"
+         solvarmod=1
+         call getin("solvarmod",solvarmod)
+         write(*,*) " solvarmod = ",solvarmod
+
          write(*,*) "date for solar flux calculation:",
-     &   " (1985 < date < 2002)"
+     &   " (1985 < date < 2002)",
+     $   " (Only used if solvarmod=0)"
          write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
          solarcondate=1993.4 ! default value
@@ -591 +599 @@
          write(*,*) " solarcondate = ",solarcondate
          
+         write(*,*) "Solar variability as observed for MY: "
+         write(*,*) "Only if solvarmod=1"
+         solvaryear=24
+         call getin("solvaryear",solvaryear)
+         write(*,*) " solvaryear = ",solvaryear
+
          write(*,*) "UV heating efficiency:",
      &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",

trunk/LMDZ.MARS/libf/phymars/physiq.F

@@ -308 +308 @@
       real rho(ngridmx,nlayermx)  ! density
       real vmr(ngridmx,nlayermx)  ! volume mixing ratio
+      real rhopart(ngridmx,nlayermx) ! number density of a given species
       real colden(ngridmx,nqmx)   ! vertical column of tracers
       REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
@@ -429 +430 @@
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-         if (callthermos) call param_read
+         if (callthermos) then
+            if(solvarmod.eq.0) call param_read
+            if(solvarmod.eq.1) call param_read_e107
+         endif
 #endif
 c        Initialize R and Cp as constant
@@ -561 +565 @@
                  CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte)
               else if(nltemodel.eq.2) then
-                 do ig=1,ngrid
-                    do l=1,nlayer
-                       co2vmr_gcm(ig,l)=pq(ig,l,igcm_co2)*
-     $                      mmean(ig,l)/mmol(igcm_co2)
-                       n2vmr_gcm(ig,l)=pq(ig,l,igcm_n2)*
-     $                      mmean(ig,l)/mmol(igcm_n2)
-                       covmr_gcm(ig,l)=pq(ig,l,igcm_co)*
-     $                      mmean(ig,l)/mmol(igcm_co)
-                       ovmr_gcm(ig,l)=pq(ig,l,igcm_o)*
-     $                      mmean(ig,l)/mmol(igcm_o)
-                    enddo
-                 enddo
+                co2vmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_co2)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co2)
+                n2vmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_n2)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_n2)
+                covmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_co)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_co)
+                ovmr_gcm(1:ngrid,1:nlayer)=
+     &                      pq(1:ngrid,1:nlayer,igcm_o)*
+     &                      mmean(1:ngrid,1:nlayer)/mmol(igcm_o)
                  
                  CALL NLTEdlvr09_TCOOL(ngrid,nlayer,pplay*9.869e-6,
@@ -578 +582 @@
      $                ovmr_gcm,  zdtnlte )
 
-                 do ig=1,ngrid
-                    do l=1,nlayer
-                       zdtnlte(ig,l)=zdtnlte(ig,l)/86400.
-                    enddo
-                 enddo
+                 zdtnlte(1:ngrid,1:nlayer)=
+     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
               endif
            else
@@ -1549 +1550 @@
         call wstats(ngrid,"v","Meridional (North-South) wind",
      &                "m.s-1",3,zv)
-c           call wstats(ngrid,"w","Vertical (down-up) wind",
-c     &                "m.s-1",3,pw)
+        call wstats(ngrid,"w","Vertical (down-up) wind",
+     &                "m.s-1",3,pw)
         call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho)
-c           call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
+        call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
 c          call wstats(ngrid,"q2",
 c    &                "Boundary layer eddy kinetic energy",
@@ -1624 +1625 @@
      $                 noms(iq) .ne. "ccn_mass" .and.
      $                 noms(iq) .ne. "ccn_number") then
-                   do l=1,nlayer
-                      do ig=1,ngrid
-                         vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
-                      end do
-                   end do
+                    vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
+                    rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
+     &                          *rho(1:ngrid,1:nlayer)*n_avog/
+     &                           (1000*mmol(iq))
                    call wstats(ngrid,"vmr_"//trim(noms(iq)),
      $                         "Volume mixing ratio","mol/mol",3,vmr)
+!                   call wstats(ngrid,"rho_"//trim(noms(iq)),
+!     $                     "Number density","cm-3",3,rhopart)
+!                   call writediagfi(ngrid,"rho_"//trim(noms(iq)),
+!     $                     "Number density","cm-3",3,rhopart)
                    if ((noms(iq).eq."o") .or. (noms(iq).eq."co2").or.
      $                 (noms(iq).eq."o3")) then                     
@@ -1928 +1933 @@
 c          endif ! (submicron)
          end if  ! (tracer.and.(dustbin.ne.0))
+
+
+c        ----------------------------------------------------------
+c        Thermospheric outputs
+c        ----------------------------------------------------------
+
+         if(callthermos) then 
+
+            call WRITEDIAGFI(ngridmx,"q15um","15 um cooling","K/s",
+     $           3,zdtnlte)
+            call WRITEDIAGFI(ngridmx,"quv","UV heating","K/s",
+     $           3,zdteuv)
+            call WRITEDIAGFI(ngridmx,"cond","Thermal conduction","K/s",
+     $           3,zdtconduc)
+            call WRITEDIAGFI(ngridmx,"qnir","NIR heating","K/s",
+     $           3,zdtnirco2)
+
+         endif  !(callthermos)
 
 c        ----------------------------------------------------------