Changeset 1232

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

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/dyn3d/dump2d.F (modified) (1 diff)
• libf/phymars/lect_start_archive.F (modified) (13 diffs)
• libf/phymars/newstart.F (modified) (8 diffs)

trunk/LMDZ.MARS/README

@@ -2045 +2045 @@
 - 'tauscaling' was missing from argument list when calling physdem1
   in testphys1d
+
+== 19/04/2014 == FF+EM
+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.
+

trunk/LMDZ.MARS/libf/dyn3d/dump2d.F

@@ -37 +37 @@
       IF(zllm.GT.zmin) THEN
       DO j=1,jm
-      WRITE(*,'(72i1)') (NINT(10.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
+!      WRITE(*,'(72i1)') (NINT(10.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
+      WRITE(*,'(72i1)') (NINT(9.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
       ENDDO
       ENDIF

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

@@ -45 +45 @@
 c=======================================================================
 
-c Old variables dimensions (from file)
-c------------------------------------
-      INTEGER   imold,jmold,lmold,nsoilold,nqold
-
-c et autres:
-c----------
-      integer,intent(in) :: ngrid ! number of atmospheric columns
+! 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" 
@@ -82 +98 @@
 
 !      INTEGER vardim(4)
-      REAL date
       INTEGER   memo
 !      character (len=50) :: tmpname
@@ -88 +103 @@
 c Variable histoire 
 c------------------
-      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
-      REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
-      REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
+      REAL ::qtot(iip1,jjp1,llm)
 
 c autre variables dynamique nouvelle grille
@@ -104 +117 @@
 c variable physique
 c------------------
-      REAL tsurf(ngrid) ! surface temperature
-      REAL tsoil(ngrid,nsoilmx) ! soil temperature
-      REAL co2ice(ngrid) ! CO2 ice layer
-      REAL emis(ngrid)
-      REAL q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
-      REAL tauscaling(ngrid) ! dust conversion factor
 c     REAL phisfi(ngrid)
 
       INTEGER i,j,l
-      INTEGER nid,nvarid
+      INTEGER nvarid
 c     REAL year_day,periheli,aphelie,peri_day
 c     REAL obliquit,z0,emin_turb,lmixmin
@@ -131 +138 @@
 c------------------------------------------------------
       real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
-      REAL phisold_newgrid(iip1,jjp1)
-      REAL t(iip1,jjp1,llm)
       real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
       real inertiedatS(iip1,jjp1,nsoilmx)
@@ -184 +189 @@
 !      real, dimension(:), allocatable :: oldmlayer
 
-      real surfith(iip1,jjp1) ! surface thermal inertia
 !      real surfithfi(ngrid)
       ! surface thermal inertia at old horizontal grid resolution
@@ -231 +235 @@
       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)
@@ -237 +247 @@
       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
 
@@ -244 +260 @@
       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
 
@@ -251 +273 @@
       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
@@ -298 +326 @@
 ! 1.3 Report dimensions
       
-      write(*,*) "Start_archive dimensions:"
+      write(*,*) "lect_start_archive: Start_archive dimensions:"
       write(*,*) "longitude: ",imold
       write(*,*) "latitude: ",jmold
       write(*,*) "altitude: ",lmold
-      write(*,*) "tracers: ",nqold
+      if (nqold.gt.0) then
+        write(*,*) "old tracers q*: ",nqold
+      endif
       write(*,*) "subsurface_layers: ",nsoilold
       if (depthinterpol) then
@@ -1281 +1311 @@
      &                   rlonuold,rlatvold,rlonu,rlatv)
       call scal_wind(us,vs,unat,vnat)
-      write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
+      write (*,*) 'lect_start_archive: unat ', unat (1,1,1)    ! INFO
       do l=1,llm
         do j = 1, jjp1
@@ -1290 +1320 @@
         end do
       end do 
-      write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
+      write (*,*) 'lect_start_archive: ucov ', ucov (1,1,1)  ! INFO
 c     write(48,*) 'ucov',ucov
       do l=1,llm

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