Changeset 1232 for trunk/LMDZ.MARS/libf/phymars/newstart.F

Timestamp:

Apr 19, 2014, 10:35:38 AM (11 years ago)

Author:

emillour

Message:

Mars GCM:
Add "composition" option in newstart to enable global modification of
CO2,N2,Ar,O2 and CO mixing ratios. Plus some cleanup on physics variables
(saved arrays) which have been moved to modules.
FF+EM

File:

• trunk/LMDZ.MARS/libf/phymars/newstart.F (modified) (8 diffs)

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

@@ -15 +15 @@
 c=======================================================================
 
-! to use  'getin'
       use ioipsl_getincom, only: getin 
       use infotrac, only: iniadvtrac, nqtot, tname
-      use tracer_mod, only: noms, igcm_dust_number, igcm_dust_mass,
+      use tracer_mod, only: noms, mmol,
+     &                      igcm_dust_number, igcm_dust_mass,
      &                      igcm_ccn_number, igcm_ccn_mass,
-     &                      igcm_h2o_vap, igcm_h2o_ice
+     &                      igcm_h2o_vap, igcm_h2o_ice, igcm_co2,
+     &                      igcm_n2, igcm_ar, igcm_o2, igcm_co
       use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
-     &                     albedodat, z0_default
-      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx
+     &                     albedodat, z0_default, qsurf, tsurf,
+     &                     co2ice, emis
+      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
       use control_mod, only: day_step, iphysiq, anneeref
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
       use comgeomphy, only: initcomgeomphy
-      use planete_h
+!      use planete_h
+      use dimradmars_mod, only: tauscaling
+      use turb_mod, only: q2
       implicit none
 
@@ -69 +73 @@
 c--------------------------------------------------
       INTEGER nid_dyn, nid_fi,nid,nvarid
-!      INTEGER size
-      INTEGER length
-      parameter (length = 100)
       INTEGER tab0
-      INTEGER   NB_ETATMAX
-      parameter (NB_ETATMAX = 100)
 
       REAL  date
@@ -105 +104 @@
 c variable physique
 c------------------
-      REAL tsurf(ngridmx)       ! surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
-      REAL co2ice(ngridmx)      ! CO2 ice layer
-      REAL emis(ngridmx)        ! surface emissivity
-      REAL tauscaling(ngridmx) ! dust conversion factor
+!      REAL tsurf(ngridmx)      ! surface temperature 
+!      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
+!      REAL co2ice(ngridmx)     ! CO2 ice layer
+!      REAL emis(ngridmx)       ! surface emissivity
+!      REAL tauscaling(ngridmx) ! dust conversion factor
       REAL tauscadyn(iip1,jjp1) ! dust conversion factor on the dynamics grid
-      REAL,ALLOCATABLE :: qsurf(:,:)
-      REAL q2(ngridmx,nlayermx+1)
+!      REAL,ALLOCATABLE :: qsurf(:,:)
+!      REAL q2(ngridmx,nlayermx+1)
 !      REAL rnaturfi(ngridmx)
       real alb(iip1,jjp1),albfi(ngridmx) ! albedos
@@ -122 +121 @@
       real mugaz ! molar mass of the atmosphere
 
-      EXTERNAL iniconst,geopot,inigeom
       integer ierr  !, nbetat
-      integer ismin
-      external ismin
 
 c Variables on the new grid along scalar points 
@@ -181 +177 @@
       real :: MONS_coeffN ! coeff for northern hemisphere
 !      real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
+! Reference position for composition
+      real :: latref,lonref,dlatmin,dlonmin
+! Variable used to change composition
+      real :: Svmr,Smmr,Smmr_old,Smmr_new,Sn
+      real :: Mair_old,Mair_new,vmr_old,vmr_new
+      real,allocatable :: coefvmr(:)  ! Correction coefficient when changing composition
 
 c sortie visu pour les champs dynamiques
@@ -200 +202 @@
 ! allocate arrays
       allocate(q(iip1,jjp1,llm,nqtot))
-      allocate(qsurf(ngridmx,nqtot))
-      
+!      allocate(qsurf(ngridmx,nqtot)) ! done in ini_surfdat_h
+      allocate(coefvmr(nqtot))
+
 c=======================================================================
 c   Choice of the start file(s) to use
@@ -521 +524 @@
       write(*,*) 'ini_q-h2o    : tracers initialization for chemistry
      $ only'
+      write(*,*) 'composition  : change atm main composition: CO2,N2,Ar,
+     $ O2,CO'
       write(*,*) 'ini_h2osurf  : reinitialize surface water ice '
       write(*,*) 'noglacier    : Remove tropical H2O ice if |lat|<45'
@@ -974 +979 @@
           write(*,*) 'inichim_newstart: chemical species initialised
      $ (except water vapour)'
+
+c      composition : change main composition: CO2,N2,Ar,O2,CO (FF 03/2014)
+c      --------------------------------------------------------
+       else if (trim(modif) .eq. 'composition') then
+          write(*,*) "Lat (degN)  lon (degE) of the reference site ?"
+          write(*,*) "e.g. MSL : -4.5  137.  "
+ 301      read(*,*,iostat=ierr) latref, lonref
+          if(ierr.ne.0) goto 301
+
+
+        !  Select GCM point close to reference site
+          dlonmin =90.
+          DO i=1,iip1-1
+             if (abs(rlonv(i)*180./pi -lonref).lt.dlonmin)then
+                iref=i
+                dlonmin=abs(rlonv(i)*180./pi -lonref)
+             end if   
+          ENDDO
+          dlatmin =45.
+          DO j=1,jjp1
+             if (abs(rlatu(j)*180./pi -latref).lt.dlatmin)then
+                jref=j
+                dlatmin=abs(rlatu(j)*180./pi -latref)
+             end if   
+          ENDDO
+          write(*,*) "In GCM : lat= " ,  rlatu(jref)*180./pi
+          write(*,*) "In GCM : lon= " ,  rlonv(iref)*180./pi
+          write(*,*)
+
+        ! Compute air molar mass at reference site
+          Smmr=0
+          Sn = 0
+          do iq=1,nqtot 
+             if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
+     &      .or. (iq.eq.igcm_ar).or.(iq.eq.igcm_o2)
+     &      .or. (iq.eq.igcm_co)) then
+                 Smmr=Smmr+q(iref,jref,1,iq)
+                 Sn=Sn+q(iref,jref,1,iq)/mmol(iq) 
+             end if
+          end do
+          write(*,*) "At reference site :  "
+          write(*,*) "Sum of mass mix. ratio (should be about 1)=",Smmr
+          Mair_old =Smmr/Sn
+          write(*,*)
+     &     "Air molar mass (g/mol) at reference site= ",Mair_old
+
+        ! Ask for new volume mixing ratio at reference site
+          Svmr =0.
+          Sn =0.
+          do iq=1,nqtot 
+           coefvmr(iq) = 1.
+           if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
+     &     .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
+
+             vmr_old=q(iref,jref,1,iq)*Mair_old/mmol(iq)  
+             write(*,*) "Previous vmr("//trim(tname(iq))//")= ", vmr_old
+
+              if (iq.eq.igcm_n2) then
+                write(*,*) "New vmr(n2)? (MSL: 1.89E-02 at Ls~186)"
+              endif
+              if (iq.eq.igcm_ar) then
+                write(*,*) "New vmr(ar)? (MSL: 1.93E-02 at Ls~186)"
+              endif
+              if (iq.eq.igcm_o2) then
+                write(*,*) "New vmr(o2)? (MSL: 1.46E-03) at Ls~186"
+              endif
+              if (iq.eq.igcm_co) then
+                write(*,*) "New vmr(co)? (~8.E-04)"
+              endif
+ 302          read(*,*,iostat=ierr) vmr_new
+              if(ierr.ne.0) goto 302
+              write(*,*) "New vmr("//trim(tname(iq))//")= ",vmr_new
+              write(*,*) 
+              coefvmr(iq) = vmr_new/vmr_old
+              Svmr=Svmr+vmr_new
+              Sn=Sn+vmr_new*mmol(iq)
+           end if
+          enddo ! of do iq=1,nqtot 
+      !  Estimation of new Air molar mass at reference site (assuming vmr_co2 = 1-Svmr)
+          Mair_new = Sn + (1-Svmr)*mmol(igcm_co2) 
+          write(*,*)
+     &     "NEW Air molar mass (g/mol) at reference site= ",Mair_new
+
+        ! Compute mass mixing ratio changes  
+          do iq=1,nqtot  
+            if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
+     &          .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
+             write(*,*) "Everywhere mmr("//trim(tname(iq))//
+     &        ") is multiplied by ",coefvmr(iq)*Mair_old/Mair_new
+            end if
+          end do
+
+        ! Recompute mass mixing ratios everywhere, and adjust mmr(CO2) to keep sum of mmr constant.
+          do l=1,llm
+           do j=1,jjp1
+            do i=1,iip1
+              Smmr_old = 0.
+              Smmr_new = 0.
+              do iq=1,nqtot  
+                if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
+     &         .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
+                   Smmr_old = Smmr_old + q(i,j,l,iq) ! sum of old mmr 
+                   q(i,j,l,iq)=q(i,j,l,iq)*coefvmr(iq)*Mair_old/Mair_new
+                   Smmr_new = Smmr_new + q(i,j,l,iq) ! sum of new mmr
+                end if 
+              enddo
+              q(i,j,l,igcm_co2) = q(i,j,l,igcm_co2) + Smmr_old-Smmr_new
+            enddo
+           enddo
+          enddo
+
+          write(*,*)
+     &   "mmr(CO2) is modified everywhere to keep sum of mmr constant"
+          write(*,*) 'At reference site vmr(CO2)=', 
+     &        q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
+          write(*,*) "Compared to MSL observation: vmr(CO2)= 0.9597" 
+
 
 c      wetstart : wet atmosphere with a north to south gradient