Changeset 1232 for trunk/LMDZ.MARS/libf


Ignore:
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/libf
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r38 r1232  
    3737      IF(zllm.GT.zmin) THEN
    3838      DO j=1,jm
    39       WRITE(*,'(72i1)') (NINT(10.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
     39!      WRITE(*,'(72i1)') (NINT(10.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
     40      WRITE(*,'(72i1)') (NINT(9.*(z(i,j)-zmin)/(zllm-zmin)),i=1,im)
    4041      ENDDO
    4142      ENDIF

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

    r1226 r1232  
    4545c=======================================================================
    4646
    47 c Old variables dimensions (from file)
    48 c------------------------------------
    49       INTEGER   imold,jmold,lmold,nsoilold,nqold
    50 
    51 c et autres:
    52 c----------
    53       integer,intent(in) :: ngrid ! number of atmospheric columns
     47! routine arguments
     48! -----------------
     49      integer,intent(in) :: ngrid ! number of atmosphferic columns
    5450                                  ! on new physics grid
    5551      integer,intent(in) :: nlayer ! number of atmospheric layers
    5652                                   ! on new grid
    5753      integer,intent(in) :: nqtot ! number of advected tracers
     54      REAL,INTENT(OUT) :: date
     55      REAL,INTENT(OUT) :: vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
     56      REAL,INTENT(OUT) :: h(iip1,jjp1,llm),ps(iip1,jjp1)
     57      REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
     58      REAL,INTENT(OUT) :: tsurf(ngrid) ! surface temperature
     59      REAL,INTENT(OUT) :: tsoil(ngrid,nsoilmx) ! soil temperature
     60      REAL,INTENT(OUT) :: co2ice(ngrid) ! CO2 ice layer
     61      REAL,INTENT(OUT) :: emis(ngrid)
     62      REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
     63      REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor
     64      REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
     65      REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
     66      real,intent(in) :: surfith(iip1,jjp1) ! surface thermal inertia
     67      INTEGER,INTENT(IN) :: nid ! input NetCDF file identifier
     68
     69
     70
     71c Old variables dimensions (from file)
     72c------------------------------------
     73      INTEGER   imold,jmold,lmold,nsoilold,nqold
    5874
    5975c Variables pour les lectures des fichiers "ini"
    … …  
    8298
    8399!      INTEGER vardim(4)
    84       REAL date
    85100      INTEGER   memo
    86101!      character (len=50) :: tmpname
    … …  
    88103c Variable histoire
    89104c------------------
    90       REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
    91       REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
    92       REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
     105      REAL ::qtot(iip1,jjp1,llm)
    93106
    94107c autre variables dynamique nouvelle grille
    … …  
    104117c variable physique
    105118c------------------
    106       REAL tsurf(ngrid) ! surface temperature
    107       REAL tsoil(ngrid,nsoilmx) ! soil temperature
    108       REAL co2ice(ngrid) ! CO2 ice layer
    109       REAL emis(ngrid)
    110       REAL q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
    111       REAL tauscaling(ngrid) ! dust conversion factor
    112119c     REAL phisfi(ngrid)
    113120
    114121      INTEGER i,j,l
    115       INTEGER nid,nvarid
     122      INTEGER nvarid
    116123c     REAL year_day,periheli,aphelie,peri_day
    117124c     REAL obliquit,z0,emin_turb,lmixmin
    … …  
    131138c------------------------------------------------------
    132139      real us(iip1,jjp1,llm),vs(iip1,jjp1,llm)
    133       REAL phisold_newgrid(iip1,jjp1)
    134       REAL t(iip1,jjp1,llm)
    135140      real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
    136141      real inertiedatS(iip1,jjp1,nsoilmx)
    … …  
    184189!      real, dimension(:), allocatable :: oldmlayer
    185190
    186       real surfith(iip1,jjp1) ! surface thermal inertia
    187191!      real surfithfi(ngrid)
    188192      ! surface thermal inertia at old horizontal grid resolution
    … …  
    231235      endif
    232236      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
     237      if (ierr.ne.NF_NOERR) then
     238        write(*,*) 'lect_start_archive error: cannot find Time length'
     239        stop
     240      else
     241        write(*,*) "lect_start_archive: timelen=",timelen
     242      endif
    233243
    234244      ierr= NF_INQ_DIMID(nid,"latitude",dimid)
    … …  
    237247      endif
    238248      ierr= NF_INQ_DIMLEN(nid,dimid,jmold)
     249      if (ierr.ne.NF_NOERR) then
     250        write(*,*) 'lect_start_archive error: cannot find lat length'
     251        stop
     252      else
     253        write(*,*) "lect_start_archive: jmold=",jmold
     254      endif
    239255      jmold=jmold-1
    240256
    … …  
    244260      endif
    245261      ierr= NF_INQ_DIMLEN(nid,dimid,imold)
     262      if (ierr.ne.NF_NOERR) then
     263        write(*,*) 'lect_start_archive error: cannot find lon length'
     264        stop
     265      else
     266        write(*,*) "lect_start_archive: imold=",imold
     267      endif
    246268      imold=imold-1
    247269
    … …  
    251273      endif
    252274      ierr= NF_INQ_DIMLEN(nid,dimid,lmold)
     275      if (ierr.ne.NF_NOERR) then
     276        write(*,*) 'lect_start_archive error: cannot find alt length'
     277        stop
     278      else
     279        write(*,*) "lect_start_archive: lmold=",lmold
     280      endif
    253281
    254282      nqold=0
    … …  
    298326! 1.3 Report dimensions
    299327     
    300       write(*,*) "Start_archive dimensions:"
     328      write(*,*) "lect_start_archive: Start_archive dimensions:"
    301329      write(*,*) "longitude: ",imold
    302330      write(*,*) "latitude: ",jmold
    303331      write(*,*) "altitude: ",lmold
    304       write(*,*) "tracers: ",nqold
     332      if (nqold.gt.0) then
     333        write(*,*) "old tracers q*: ",nqold
     334      endif
    305335      write(*,*) "subsurface_layers: ",nsoilold
    306336      if (depthinterpol) then
    … …  
    12811311     &                   rlonuold,rlatvold,rlonu,rlatv)
    12821312      call scal_wind(us,vs,unat,vnat)
    1283       write (*,*) 'lect_start_archive: unat ', unat (1,2,1)    ! INFO
     1313      write (*,*) 'lect_start_archive: unat ', unat (1,1,1)    ! INFO
    12841314      do l=1,llm
    12851315        do j = 1, jjp1
    … …  
    12901320        end do
    12911321      end do
    1292       write (*,*) 'lect_start_archive: ucov ', ucov (1,2,1)  ! INFO
     1322      write (*,*) 'lect_start_archive: ucov ', ucov (1,1,1)  ! INFO
    12931323c     write(48,*) 'ucov',ucov
    12941324      do l=1,llm

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

    r1231 r1232  
    1515c=======================================================================
    1616
    17 ! to use  'getin'
    1817      use ioipsl_getincom, only: getin
    1918      use infotrac, only: iniadvtrac, nqtot, tname
    20       use tracer_mod, only: noms, igcm_dust_number, igcm_dust_mass,
     19      use tracer_mod, only: noms, mmol,
     20     &                      igcm_dust_number, igcm_dust_mass,
    2121     &                      igcm_ccn_number, igcm_ccn_mass,
    22      &                      igcm_h2o_vap, igcm_h2o_ice
     22     &                      igcm_h2o_vap, igcm_h2o_ice, igcm_co2,
     23     &                      igcm_n2, igcm_ar, igcm_o2, igcm_co
    2324      use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
    24      &                     albedodat, z0_default
    25       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx
     25     &                     albedodat, z0_default, qsurf, tsurf,
     26     &                     co2ice, emis
     27      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
    2628      use control_mod, only: day_step, iphysiq, anneeref
    2729      use phyredem, only: physdem0, physdem1
    2830      use iostart, only: open_startphy
    2931      use comgeomphy, only: initcomgeomphy
    30       use planete_h
     32!      use planete_h
     33      use dimradmars_mod, only: tauscaling
     34      use turb_mod, only: q2
    3135      implicit none
    3236
    … …  
    6973c--------------------------------------------------
    7074      INTEGER nid_dyn, nid_fi,nid,nvarid
    71 !      INTEGER size
    72       INTEGER length
    73       parameter (length = 100)
    7475      INTEGER tab0
    75       INTEGER   NB_ETATMAX
    76       parameter (NB_ETATMAX = 100)
    7776
    7877      REAL  date
    … …  
    105104c variable physique
    106105c------------------
    107       REAL tsurf(ngridmx)       ! surface temperature
    108       REAL tsoil(ngridmx,nsoilmx) ! soil temperature
    109       REAL co2ice(ngridmx)      ! CO2 ice layer
    110       REAL emis(ngridmx)        ! surface emissivity
    111       REAL tauscaling(ngridmx) ! dust conversion factor
     106!      REAL tsurf(ngridmx)      ! surface temperature
     107!      REAL tsoil(ngridmx,nsoilmx) ! soil temperature
     108!      REAL co2ice(ngridmx)     ! CO2 ice layer
     109!      REAL emis(ngridmx)       ! surface emissivity
     110!      REAL tauscaling(ngridmx) ! dust conversion factor
    112111      REAL tauscadyn(iip1,jjp1) ! dust conversion factor on the dynamics grid
    113       REAL,ALLOCATABLE :: qsurf(:,:)
    114       REAL q2(ngridmx,nlayermx+1)
     112!      REAL,ALLOCATABLE :: qsurf(:,:)
     113!      REAL q2(ngridmx,nlayermx+1)
    115114!      REAL rnaturfi(ngridmx)
    116115      real alb(iip1,jjp1),albfi(ngridmx) ! albedos
    … …  
    122121      real mugaz ! molar mass of the atmosphere
    123122
    124       EXTERNAL iniconst,geopot,inigeom
    125123      integer ierr  !, nbetat
    126       integer ismin
    127       external ismin
    128124
    129125c Variables on the new grid along scalar points
    … …  
    181177      real :: MONS_coeffN ! coeff for northern hemisphere
    182178!      real,parameter :: icedepthmin=1.e-3 ! Ice begins at most at that depth
     179! Reference position for composition
     180      real :: latref,lonref,dlatmin,dlonmin
     181! Variable used to change composition
     182      real :: Svmr,Smmr,Smmr_old,Smmr_new,Sn
     183      real :: Mair_old,Mair_new,vmr_old,vmr_new
     184      real,allocatable :: coefvmr(:)  ! Correction coefficient when changing composition
    183185
    184186c sortie visu pour les champs dynamiques
    … …  
    200202! allocate arrays
    201203      allocate(q(iip1,jjp1,llm,nqtot))
    202       allocate(qsurf(ngridmx,nqtot))
    203      
     204!      allocate(qsurf(ngridmx,nqtot)) ! done in ini_surfdat_h
     205      allocate(coefvmr(nqtot))
     206
    204207c=======================================================================
    205208c   Choice of the start file(s) to use
    … …  
    521524      write(*,*) 'ini_q-h2o    : tracers initialization for chemistry
    522525     $ only'
     526      write(*,*) 'composition  : change atm main composition: CO2,N2,Ar,
     527     $ O2,CO'
    523528      write(*,*) 'ini_h2osurf  : reinitialize surface water ice '
    524529      write(*,*) 'noglacier    : Remove tropical H2O ice if |lat|<45'
    … …  
    974979          write(*,*) 'inichim_newstart: chemical species initialised
    975980     $ (except water vapour)'
     981
     982c      composition : change main composition: CO2,N2,Ar,O2,CO (FF 03/2014)
     983c      --------------------------------------------------------
     984       else if (trim(modif) .eq. 'composition') then
     985          write(*,*) "Lat (degN)  lon (degE) of the reference site ?"
     986          write(*,*) "e.g. MSL : -4.5  137.  "
     987 301      read(*,*,iostat=ierr) latref, lonref
     988          if(ierr.ne.0) goto 301
     989
     990
     991        !  Select GCM point close to reference site
     992          dlonmin =90.
     993          DO i=1,iip1-1
     994             if (abs(rlonv(i)*180./pi -lonref).lt.dlonmin)then
     995                iref=i
     996                dlonmin=abs(rlonv(i)*180./pi -lonref)
     997             end if   
     998          ENDDO
     999          dlatmin =45.
     1000          DO j=1,jjp1
     1001             if (abs(rlatu(j)*180./pi -latref).lt.dlatmin)then
     1002                jref=j
     1003                dlatmin=abs(rlatu(j)*180./pi -latref)
     1004             end if   
     1005          ENDDO
     1006          write(*,*) "In GCM : lat= " ,  rlatu(jref)*180./pi
     1007          write(*,*) "In GCM : lon= " ,  rlonv(iref)*180./pi
     1008          write(*,*)
     1009
     1010        ! Compute air molar mass at reference site
     1011          Smmr=0
     1012          Sn = 0
     1013          do iq=1,nqtot
     1014             if ((iq.eq.igcm_co2).or.(iq.eq.igcm_n2)
     1015     &      .or. (iq.eq.igcm_ar).or.(iq.eq.igcm_o2)
     1016     &      .or. (iq.eq.igcm_co)) then
     1017                 Smmr=Smmr+q(iref,jref,1,iq)
     1018                 Sn=Sn+q(iref,jref,1,iq)/mmol(iq)
     1019             end if
     1020          end do
     1021          write(*,*) "At reference site :  "
     1022          write(*,*) "Sum of mass mix. ratio (should be about 1)=",Smmr
     1023          Mair_old =Smmr/Sn
     1024          write(*,*)
     1025     &     "Air molar mass (g/mol) at reference site= ",Mair_old
     1026
     1027        ! Ask for new volume mixing ratio at reference site
     1028          Svmr =0.
     1029          Sn =0.
     1030          do iq=1,nqtot
     1031           coefvmr(iq) = 1.
     1032           if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
     1033     &     .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
     1034
     1035             vmr_old=q(iref,jref,1,iq)*Mair_old/mmol(iq) 
     1036             write(*,*) "Previous vmr("//trim(tname(iq))//")= ", vmr_old
     1037
     1038              if (iq.eq.igcm_n2) then
     1039                write(*,*) "New vmr(n2)? (MSL: 1.89E-02 at Ls~186)"
     1040              endif
     1041              if (iq.eq.igcm_ar) then
     1042                write(*,*) "New vmr(ar)? (MSL: 1.93E-02 at Ls~186)"
     1043              endif
     1044              if (iq.eq.igcm_o2) then
     1045                write(*,*) "New vmr(o2)? (MSL: 1.46E-03) at Ls~186"
     1046              endif
     1047              if (iq.eq.igcm_co) then
     1048                write(*,*) "New vmr(co)? (~8.E-04)"
     1049              endif
     1050 302          read(*,*,iostat=ierr) vmr_new
     1051              if(ierr.ne.0) goto 302
     1052              write(*,*) "New vmr("//trim(tname(iq))//")= ",vmr_new
     1053              write(*,*)
     1054              coefvmr(iq) = vmr_new/vmr_old
     1055              Svmr=Svmr+vmr_new
     1056              Sn=Sn+vmr_new*mmol(iq)
     1057           end if
     1058          enddo ! of do iq=1,nqtot
     1059      !  Estimation of new Air molar mass at reference site (assuming vmr_co2 = 1-Svmr)
     1060          Mair_new = Sn + (1-Svmr)*mmol(igcm_co2)
     1061          write(*,*)
     1062     &     "NEW Air molar mass (g/mol) at reference site= ",Mair_new
     1063
     1064        ! Compute mass mixing ratio changes 
     1065          do iq=1,nqtot 
     1066            if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
     1067     &          .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
     1068             write(*,*) "Everywhere mmr("//trim(tname(iq))//
     1069     &        ") is multiplied by ",coefvmr(iq)*Mair_old/Mair_new
     1070            end if
     1071          end do
     1072
     1073        ! Recompute mass mixing ratios everywhere, and adjust mmr(CO2) to keep sum of mmr constant.
     1074          do l=1,llm
     1075           do j=1,jjp1
     1076            do i=1,iip1
     1077              Smmr_old = 0.
     1078              Smmr_new = 0.
     1079              do iq=1,nqtot 
     1080                if ((iq.eq.igcm_n2).or.(iq.eq.igcm_ar)
     1081     &         .or. (iq.eq.igcm_o2).or.(iq.eq.igcm_co)) then
     1082                   Smmr_old = Smmr_old + q(i,j,l,iq) ! sum of old mmr
     1083                   q(i,j,l,iq)=q(i,j,l,iq)*coefvmr(iq)*Mair_old/Mair_new
     1084                   Smmr_new = Smmr_new + q(i,j,l,iq) ! sum of new mmr
     1085                end if
     1086              enddo
     1087              q(i,j,l,igcm_co2) = q(i,j,l,igcm_co2) + Smmr_old-Smmr_new
     1088            enddo
     1089           enddo
     1090          enddo
     1091
     1092          write(*,*)
     1093     &   "mmr(CO2) is modified everywhere to keep sum of mmr constant"
     1094          write(*,*) 'At reference site vmr(CO2)=',
     1095     &        q(iref,jref,1,igcm_co2)*Mair_new/mmol(igcm_co2)
     1096          write(*,*) "Compared to MSL observation: vmr(CO2)= 0.9597"
     1097
    9761098
    9771099c      wetstart : wet atmosphere with a north to south gradient
Note: See TracChangeset for help on using the changeset viewer.