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

Timestamp:

Mar 2, 2012, 9:53:27 AM (13 years ago)

Author:

emillour

Message:

Mars GCM:

bug fix in "aeronomars/moldiff_red.F90" (wrong mmol array size in included

subroutine) + adapted extreme limit test for H to altitudes of ~4000km
(compatble with 50 layer model).

bug fix in "nirco2abs.F" index of "co2" and "o" tracers were hard coded as

1 and 3!

updates from FGG of euvheat.F, callkeys.h and inifis.F to have the

"euveff" parameter read from callphys.def

EM

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• callkeys.h (modified) (2 diffs)
• inifis.F (modified) (2 diffs)
• nirco2abs.F (modified) (6 diffs)
• tabfi.F (modified) (18 diffs)

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

@@ -18 +18 @@
      &   ,dustbin,nqchem_min,nltemodel,nircorr
      
-      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,            &
+      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan,euveff,     &
      &   tke_heat_flux
      
@@ -37 +37 @@
       real alphan
       real solarcondate
+      real euveff
       real tke_heat_flux
 

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

@@ -249 +249 @@
      $              "1-> new correction",
      $              "(matters only if callnirco2=T)"
-         nircorr=0
+         nircorr=1      !default value
          call getin("nircorr",nircorr)
          write(*,*) " nircorr = ",nircorr
@@ -594 +594 @@
          write(*,*) " solarcondate = ",solarcondate
          
+         write(*,*) "UV heating efficiency:",
+     &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
+     &   "lower values may be used to compensate low 15 um cooling"
+         euveff=0.21 !default value
+         call getin("euveff",euveff)
+         write(*,*) " euveff = ", euveff
 
          if (.not.callthermos) then

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

@@ -50 +50 @@
 #include "comdiurn.h"
 #include "nirdata.h"
+#include "tracer.h"
 
 c-----------------------------------------------------------------------
 c    Input/Output
 c    ------------
-      INTEGER ngrid,nlayer
-
-      REAL pplay(ngrid,nlayer)
-      REAL dist_sol,mu0(ngridmx),fract(ngridmx),declin
-
-      REAL pdtnirco2(ngrid,nlayer)
+      integer,intent(in) :: ngrid ! number of (horizontal) grid points
+      integer,intent(in) :: nlayer ! number of atmospheric layers
+      real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
+      real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
+      integer,intent(in) :: nq ! number of tracers
+      real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
+      real,intent(in) :: mu0(ngridmx) ! solar angle
+      real,intent(in) :: fract(ngridmx) ! day fraction of the time interval
+      real,intent(in) :: declin ! latitude of sub-solar point
+      
+      real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
 c
 c    Local variables :
@@ -73 +79 @@
 c   local saved variables
 c   ---------------------
-
+      logical,save :: firstcall=.true.
+      real,save :: ico2=0 ! index of "co2" tracer
+      real,save :: io=0 ! index of "o" tracer
 c     p0noonlte is a pressure below which non LTE effects are significant.
 c     REAL p0nonlte
@@ -90 +98 @@
       real    p2011,cociente1,merge
       real    cor0,oco2gcm
-      integer nq
-      real    pq(ngrid,nlayer,nq)
 
 c----------------------------------------------------------------------
@@ -97 +103 @@
 c     Initialisation
 c     --------------
+      if (firstcall) then
+        if (nircorr.eq.1) then
+          ! we will need co2 and o tracers
+          ico2=igcm_co2
+          if (ico2==0) then
+            write(*,*) "nirco2abs error: I need a CO2 tracer"
+            write(*,*) "     when running with nircorr==1"
+            stop
+          endif
+          io=igcm_o
+          if (io==0) then
+            write(*,*) "nirco2abs error: I need an O tracer"
+            write(*,*) "     when running with nircorr==1"
+            stop
+          endif
+        endif
+        firstcall=.false.
+      endif
+
+
 c     co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
       co2heat0=n_a*(1.52/dist_sol)**2/daysec
@@ -121 +147 @@
                if(nircorr.eq.1) then
                   cor0=1./(1.+n_p0/pplay(ig,l))**n_b
-                  if(pq(ig,l,1).gt.1.e-6) then
-                     oco2gcm=pq(ig,l,3)/pq(ig,l,1)
+                  if(pq(ig,l,ico2).gt.1.e-6) then
+                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                   else
                      oco2gcm=1.e6
@@ -179 +205 @@
                   if(nircorr.eq.1) then
                      cor0=1./(1.+n_p0/pplay(ig,l))**n_b
-                     oco2gcm=pq(ig,l,3)/pq(ig,l,1)
+                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
                      cociente1=oco2gcm/oldoco2(l)
                      merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*

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

@@ -81 +81 @@
 c -----------------------
       INTEGER setname, cluvdb, getdat
-
-      INTEGER lnblnk
 
 c-----------------------------------------------------------------------
@@ -296 +294 @@
  
  
-      do while(modif(1:1).ne.'hello')
+      do ! neverending loop
         write(*,*)
         write(*,*)
@@ -306 +304 @@
  
         write(*,*)
-        write(*,*) modif(1:lnblnk(modif)) , ' : '
-
-        if (modif(1:lnblnk(modif)) .eq. 'day_ini') then
+        write(*,*) trim(modif) , ' : '
+
+        if (trim(modif) .eq. 'day_ini') then
           write(*,*) 'current value:',day_ini
           write(*,*) 'enter new value:'
@@ -316 +314 @@
           write(*,*) 'day_ini (new value):',day_ini
 
-        else if (modif(1:lnblnk(modif)) .eq. 'z0') then
+        else if (trim(modif) .eq. 'z0') then
           write(*,*) 'current value (m):',z0_default
           write(*,*) 'enter new value (m):'
@@ -324 +322 @@
           write(*,*) ' z0 (new value):',z0_default
 
-        else if (modif(1:lnblnk(modif)) .eq. 'emin_turb') then
+        else if (trim(modif) .eq. 'emin_turb') then
           write(*,*) 'current value:',emin_turb
           write(*,*) 'enter new value:'
@@ -332 +330 @@
           write(*,*) ' emin_turb (new value):',emin_turb
 
-        else if (modif(1:lnblnk(modif)) .eq. 'lmixmin') then
+        else if (trim(modif) .eq. 'lmixmin') then
           write(*,*) 'current value:',lmixmin
           write(*,*) 'enter new value:'
@@ -340 +338 @@
           write(*,*) ' lmixmin (new value):',lmixmin
 
-        else if (modif(1:lnblnk(modif)) .eq. 'emissiv') then
+        else if (trim(modif) .eq. 'emissiv') then
           write(*,*) 'current value:',emissiv
           write(*,*) 'enter new value:'
@@ -348 +346 @@
           write(*,*) ' emissiv (new value):',emissiv
 
-        else if (modif(1:lnblnk(modif)) .eq. 'emisice') then
+        else if (trim(modif) .eq. 'emisice') then
           write(*,*) 'current value emisice(1) North:',emisice(1)
           write(*,*) 'enter new value:'
@@ -364 +362 @@
           write(*,*) ' emisice(2) (new value):',emisice(2)
 
-        else if (modif(1:lnblnk(modif)) .eq. 'albedice') then
+        else if (trim(modif) .eq. 'albedice') then
           write(*,*) 'current value albedice(1) North:',albedice(1)
           write(*,*) 'enter new value:'
@@ -380 +378 @@
           write(*,*) ' albedice(2) (new value):',albedice(2)
 
-        else if (modif(1:lnblnk(modif)) .eq. 'iceradius') then
+        else if (trim(modif) .eq. 'iceradius') then
           write(*,*) 'current value iceradius(1) North:',iceradius(1)
           write(*,*) 'enter new value:'
@@ -396 +394 @@
           write(*,*) ' iceradius(2) (new value):',iceradius(2)
 
-        else if (modif(1:lnblnk(modif)) .eq. 'dtemisice') then
+        else if (trim(modif) .eq. 'dtemisice') then
           write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
           write(*,*) 'enter new value:'
@@ -412 +410 @@
           write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
 
-        else if (modif(1:lnblnk(modif)) .eq. 'tauvis') then
+        else if (trim(modif) .eq. 'tauvis') then
           write(*,*) 'current value:',tauvis
           write(*,*) 'enter new value:'
@@ -420 +418 @@
           write(*,*) ' tauvis (new value):',tauvis
 
-        else if (modif(1:lnblnk(modif)) .eq. 'obliquit') then
+        else if (trim(modif) .eq. 'obliquit') then
           write(*,*) 'current value:',obliquit
           write(*,*) 'obliquit should be 25.19 on current Mars'
@@ -429 +427 @@
           write(*,*) ' obliquit (new value):',obliquit
 
-        else if (modif(1:lnblnk(modif)) .eq. 'peri_day') then
+        else if (trim(modif) .eq. 'peri_day') then
           write(*,*) 'current value:',peri_day
           write(*,*) 'peri_day should be 485 on current Mars'
@@ -438 +436 @@
           write(*,*) ' peri_day (new value):',peri_day
 
-        else if (modif(1:lnblnk(modif)) .eq. 'periheli') then
+        else if (trim(modif) .eq. 'periheli') then
           write(*,*) 'current value:',periheli
           write(*,*) 'perihelion should be 206.66 on current Mars'
@@ -447 +445 @@
           write(*,*) ' periheli (new value):',periheli
  
-        else if (modif(1:lnblnk(modif)) .eq. 'aphelie') then
+        else if (trim(modif) .eq. 'aphelie') then
           write(*,*) 'current value:',aphelie
           write(*,*) 'aphelion should be 249.22 on current Mars'
@@ -456 +454 @@
           write(*,*) ' aphelie (new value):',aphelie
  
-        else if (modif(1:lnblnk(modif)) .eq. 'volcapa') then
+        else if (trim(modif) .eq. 'volcapa') then
           write(*,*) 'current value:',volcapa
           write(*,*) 'enter new value:'
@@ -465 +463 @@
  
         endif
-      enddo ! of do while(modif(1:1).ne.'hello')
+      enddo ! of do ! neverending loop
 
  999  continue