Changeset 589 for trunk


Ignore:
Timestamp:
Mar 19, 2012, 4:44:04 PM (13 years ago)
Author:
aslmd
Message:

LMDZ.GENERIC: Cleaned rcm1d.F and made it truly generic by asking for planetary constants without any default values. Settings are now needed in a rcm1d.def file. Unbeknown to the user, we create a minimal run.def file, read parameters, then remove this dummy run.def. Introduced a keyword force_cpp if the user wants to give values for cpp and mugaz in def files.

Location:
trunk/LMDZ.GENERIC
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r588 r589  
    610610  in interp_horiz (which "fixes" an odd behaviour which fills some arrays with
    611611  zeros, but only when using some versions of ifort!)
     612
     613== 19/03/2012 == AS
     614- Cleaned rcm1d.F and made it truly generic by asking for planetary constants without any default values.
     615- Settings are now needed in a rcm1d.def file. Unbeknown to the user, we create a minimal run.def file, read parameters, then remove this dummy run.def.
     616- Introduced a keyword force_cpp if the user wants to give values for cpp and mugaz in def files.

  • trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90

    r588 r589  
    99!     unless option 'check_cpp_match' is set to false in
    1010!     callphys.def.
    11 !     NOTE: for now, in 1D we do as before. Jeremy, if you're
    12 !     re-writing rcm1d you may want to alter this.
    1311!
    1412!     Authors
    1513!     -------
    1614!     Robin Wordsworth (2009)
     15!     A. Spiga: make the routine OK with latest changes in rcm1d
    1716!
    1817!==================================================================
    … …  
    2625#include "callkeys.h"
    2726
    28       real cpp_c
     27      real cpp_c   
    2928      real mugaz_c
    3029
    … …  
    7473      print*,'Predefined Mg in physics is ',mugaz,'amu'
    7574
    76       if(ngridmx.eq.1)then
    77          print*,'Automatically setting cpp & mugaz to calculated values in calc_cpp_mugaz'
    78          cpp   = cpp_c
    79          mugaz = mugaz_c
    80          R     = 8.314511E+0 *1000.E+0/mugaz
    81          rcp   = R/cpp
    82 !      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
    83          ! Ehouarn: tolerate a small mismatch between computed/stored values
    84        elseif((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
     75      if (check_cpp_match) then
     76         print*,'REQUEST TO CHECK cpp_match :'
     77         if((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
    8578              (abs(1.-mugaz/mugaz_c).gt.1.e-6)) then
    86          if(check_cpp_match)then
    87             print*,'Values do not match!'
    88             print*,'Either adjust cpp / mugaz via newstart to calculated values,'
    89             print*,'or set check_cpp_match to .false. in callphys.def.'
     79            ! Ehouarn: tolerate a small mismatch between computed/stored values
     80            print*,'--> Values do not match!'
     81            print*,'    Either adjust cpp / mugaz via newstart to calculated values,'
     82            print*,'    or set check_cpp_match to .false. in callphys.def.'
    9083            stop
     84         else
     85            print*,'--> OK. Settings match composition.'
    9186         endif
    9287      endif
    9388
     89      if (.not.force_cpp) then
     90          print*,'*** Setting cpp & mugaz to computations in calc_cpp_mugaz.'
     91          mugaz = mugaz_c
     92          cpp = cpp_c
     93      else
     94          print*,'*** Setting cpp & mugaz to predefined values.'
     95      endif
     96
     97
    9498      return
    9599    end subroutine calc_cpp_mugaz

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

    r538 r589  
    2121     &   , newtonian                                                    &
    2222     &   , check_cpp_match                                              &
     23     &   , force_cpp                                                    &
    2324     &   , tau_relax                                                    &
    2425     &   , testradtimes                                                 &
    … …  
    6364      logical newtonian
    6465      logical check_cpp_match
     66      logical force_cpp
    6567      logical testradtimes
    6668      logical rayleigh

  • trunk/LMDZ.GENERIC/libf/phystd/inifis.F

    r586 r589  
    182182         write(*,*) " check_cpp_match = ",check_cpp_match
    183183
    184 
    185184         write(*,*) "call radiative transfer ?"
    186185         callrad=.true. ! default value
    … …  
    308307           call abort
    309308         endif
    310          if (noradsurf.and.calldifv) then
    311            print*,'noradsurf not compatible with a boundary layer!'
    312            call abort
    313          endif
     309         !if (noradsurf.and.calldifv) then
     310         !  print*,'noradsurf not compatible with a boundary layer!'
     311         !  call abort
     312         !endif
    314313
    315314         write(*,*)"Use Newtonian cooling for radiative transfer?"
    … …  
    508507         endif
    509508
    510          mugaz=8.314*1000./pr
     509         write(*,*) "Does user want to force cpp and mugaz?"
     510         force_cpp=.false. ! default value
     511         call getin("force_cpp",force_cpp)
     512         write(*,*) " force_cpp = ",force_cpp
     513
     514         if (force_cpp) then
     515           mugaz = -99999.
     516           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
     517           call getin("mugaz",mugaz)
     518           IF (mugaz.eq.-99999.) THEN
     519               PRINT *, "mugaz must be set if force_cpp = T"
     520               STOP
     521           ELSE
     522               write(*,*) "mugaz=",mugaz
     523           ENDIF
     524           cpp = -99999.
     525           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
     526           call getin("cpp",cpp)
     527           IF (cpp.eq.-99999.) THEN
     528               PRINT *, "cpp must be set if force_cpp = T"
     529               STOP
     530           ELSE
     531               write(*,*) "cpp=",cpp
     532           ENDIF
     533         else
     534           mugaz=8.314*1000./pr
     535         endif
    511536         call su_gases
    512537         call calc_cpp_mugaz

  • trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F

    r586 r589  
    2525!     F. Montmessin (water ice added)
    2626!     R. Wordsworth
     27!     B. Charnay
     28!     A. Spiga
    2729!
    2830!==================================================================
    … …  
    116118      saveprofile=.false.
    117119
    118 c ------------------------------------------------------
    119 c  Default values for constants (corresponding to Mars)
    120 c ------------------------------------------------------
     120c ----------------------------------------
     121c  Default values (corresponding to Mars)
     122c ----------------------------------------
    121123
    122124      pi=2.E+0*asin(1.E+0)
    123125
    124 c     Constante de la Planete Mars
    125 c     ----------------------------
    126       rad=3397200.               ! rayon de mars (m)  ~3397200 m
    127       daysec=88775.              ! duree du sol (s)  ~88775 s
    128       omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
    129       g=3.72                     ! gravite (m.s-2) ~3.72 
    130       mugaz=43.49                ! Masse molaire de l'atm (g.mol-1) ~43.49
    131       rcp=.256793                ! = r/cp  ~0.256793
    132       r= 8.314511E+0 *1000.E+0/mugaz
    133       cpp= r/rcp
    134       year_day = 669           ! duree de l'annee (sols) ~668.6
    135       periastr = 206.66          ! min. dist. star-planet (Mkm) ~206.66
    136       apoastr = 249.22           ! max. dist. star-planet (Mkm) ~249.22
    137       peri_day =  485.           ! date du periastron (sols depuis printemps)
    138       obliquit = 25.2           ! Obliquite de la planete (deg) ~25.2         
    139  
    140126c     Parametres Couche limite et Turbulence
    141127c     --------------------------------------
    … …  
    159145
    160146c ------------------------------------------------------
    161 c  Load parameters from "run.def"
     147c  Load parameters from "run.def" and "gases.def"
    162148c ------------------------------------------------------
    163149
    164 ! also check if 'run1d.def' file is around (otherwise reading parameters
    165 ! from callphys.def via getin() routine won't work.
    166       open(90,file='run.def',status='old',form='formatted',
     150! check if 'rcm1d.def' file is around
     151      open(90,file='rcm1d.def',status='old',form='formatted',
    167152     &     iostat=ierr)
    168153      if (ierr.ne.0) then
    169         write(*,*) 'Cannot find required file "run.def"'
    170         write(*,*) '  (which should contain some input parameters'
    171         write(*,*) '   along with the following line:'
    172         write(*,*) '   INCLUDEDEF=callphys.def'
    173         write(*,*) '   )'
     154        write(*,*) 'Cannot find required file "rcm1d.def"'
     155        write(*,*) 'which should contain some input parameters'
    174156        write(*,*) ' ... might as well stop here ...'
    175157        stop
    … …  
    178160      endif
    179161
     162! now, run.def is needed anyway. so we create a dummy temporary one
     163! for ioipsl to work. if a run.def is already here, stop the
     164! program and ask the user to do a bit of cleaning
     165      open(90,file='run.def',status='old',form='formatted',
     166     &     iostat=ierr)
     167      if (ierr.eq.0) then
     168        close(90)
     169        write(*,*) 'There is already a run.def file.'
     170        write(*,*) 'This is not compatible with 1D runs.'
     171        write(*,*) 'Please remove the file and restart the run.'
     172        write(*,*) 'Runtime parameters are supposed to be in rcm1d.def'
     173        stop
     174      else
     175        call system('touch run.def')
     176        call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
     177        call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
     178      endif
     179
    180180! check if we are going to run with or without tracers
    181181      write(*,*) "Run with or without tracer transport ?"
    … …  
    183183      call getin("tracer",tracer)
    184184      write(*,*) " tracer = ",tracer
     185
     186! OK. now that run.def has been read once -- any variable is in memory.
     187! so we can dump the dummy run.def
     188      call system("rm -rf run.def")
    185189
    186190! while we're at it, check if there is a 'traceur.def' file
    … …  
    256260      endif ! of if (tracer)
    257261
     262!!! We have to check that check_cpp_match is F for 1D computations
     263!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
     264      check_cpp_match = .false.
     265      call getin("check_cpp_match",check_cpp_match)
     266      if (check_cpp_match) then
     267          print*,"In 1D modeling, check_cpp_match is supposed to be F"
     268          print*,"Please correct callphys.def"
     269          stop
     270      endif
     271
     272!!! GEOGRAPHICAL INITIALIZATIONS
     273     !!! AREA. useless in 1D
     274      area(1)=1.E+0
     275      aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
     276     !!! GEOPOTENTIAL. useless in 1D because control bu surface pressure
     277      phisfi(1)=0.E+0
     278     !!! LATITUDE. can be set.
     279      lati(1)=0 ! default value for lati(1)
     280      PRINT *,'latitude (in degrees) ?'
     281      call getin("latitude",lati(1))
     282      write(*,*) " latitude = ",lati(1)
     283      lati(1)=lati(1)*pi/180.E+0
     284     !!! LONGITUDE. useless in 1D.
     285      long(1)=0.E+0
     286      long(1)=long(1)*pi/180.E+0
     287
     288!!! CALL INIFIS
     289!!! - read callphys.def
     290!!! - calculate sine and cosine of longitude and latitude
     291!!! - calculate mugaz and cp
     292!!! NB: some operations are being done dummily in inifis in 1D
     293!!! - physical constants: nevermind, things are done allright below
     294!!! - physical frequency: nevermind, in inifis this is a simple print
     295      CALL inifis(1,llm,day0,daysec,dtphys,
     296     .            lati,long,area,rad,g,r,cpp)
     297!!! We check everything is OK.
     298      PRINT *,"CHECK"
     299      PRINT *,"--> mugaz = ",mugaz
     300      PRINT *,"--> cpp = ",cpp
     301      r = 8.314511E+0 * 1000.E+0 / mugaz
     302      rcp = r / cpp
     303
     304
     305!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     306!!!! PLANETARY CONSTANTS !!!!
     307!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     308
     309      g = -99999.
     310      PRINT *,'GRAVITY in m s-2 ?'
     311      call getin("g",g)
     312      IF (g.eq.-99999.) THEN
     313          PRINT *,"STOP. I NEED g IN RUN.DEF."
     314          STOP
     315      ELSE
     316          PRINT *,"--> g = ",g
     317      ENDIF
     318
     319      !rad = -99999.
     320      !PRINT *,'PLANETARY RADIUS in m ?'
     321      !call getin("rad",rad)
     322      !IF (rad.eq.-99999.) THEN
     323      !    PRINT *,"STOP. I NEED rad IN RUN.DEF."
     324      !    STOP
     325      !ELSE
     326      !    PRINT *,"--> rad = ",rad
     327      !ENDIF
     328
     329      daysec = -99999.
     330      PRINT *,'LENGTH OF A DAY in s ?'
     331      call getin("daysec",daysec)
     332      IF (daysec.eq.-99999.) THEN
     333          PRINT *,"STOP. I NEED daysec IN RUN.DEF."
     334          STOP
     335      ELSE
     336          PRINT *,"--> daysec = ",daysec
     337      ENDIF
     338      omeg=4.*asin(1.)/(daysec)
     339      PRINT *,"OK. FROM THIS I WORKED OUT:"
     340      PRINT *,"--> omeg = ",omeg
     341
     342      year_day = -99999.
     343      PRINT *,'LENGTH OF A YEAR in days ?'
     344      call getin("year_day",year_day)
     345      IF (year_day.eq.-99999.) THEN
     346          PRINT *,"STOP. I NEED year_day IN RUN.DEF."
     347          STOP
     348      ELSE
     349          PRINT *,"--> year_day = ",year_day
     350      ENDIF
     351
     352      periastr = -99999.
     353      PRINT *,'MIN DIST STAR-PLANET in AU ?'
     354      call getin("periastr",periastr)
     355      IF (periastr.eq.-99999.) THEN
     356          PRINT *,"STOP. I NEED periastr IN RUN.DEF."
     357          STOP
     358      ELSE
     359          PRINT *,"--> periastr = ",periastr
     360      ENDIF
     361      periastr=periastr*149.598 ! AU to Mkm
     362
     363      apoastr = -99999.
     364      PRINT *,'MIN DIST STAR-PLANET in AU ?'
     365      call getin("apoastr",apoastr)
     366      IF (apoastr.eq.-99999.) THEN
     367          PRINT *,"STOP. I NEED apoastr IN RUN.DEF."
     368          STOP
     369      ELSE
     370          PRINT *,"--> apoastr = ",apoastr
     371      ENDIF
     372      apoastr=apoastr*149.598 ! AU to Mkm
     373
     374      peri_day = -99999.
     375      PRINT *,'DATE OF PERIASTRON in days ?'
     376      call getin("peri_day",peri_day)
     377      IF (peri_day.eq.-99999.) THEN
     378          PRINT *,"STOP. I NEED peri_day IN RUN.DEF."
     379          STOP
     380      ELSE IF (peri_day.gt.year_day) THEN
     381          PRINT *,"STOP. peri_day.gt.year_day"
     382          STOP
     383      ELSE 
     384          PRINT *,"--> peri_day = ", peri_day
     385      ENDIF
     386
     387      obliquit = -99999.
     388      PRINT *,'OBLIQUITY in deg ?'
     389      call getin("obliquit",obliquit)
     390      IF (obliquit.eq.-99999.) THEN
     391          PRINT *,"STOP. I NEED obliquit IN RUN.DEF."
     392          STOP
     393      ELSE
     394          PRINT *,"--> obliquit = ",obliquit
     395      ENDIF
     396
     397      psurf = -99999.
     398      PRINT *,'SURFACE PRESSURE in Pa ?'
     399      call getin("psurf",psurf)
     400      IF (psurf.eq.-99999.) THEN
     401          PRINT *,"STOP. I NEED psurf IN RUN.DEF."
     402          STOP
     403      ELSE
     404          PRINT *,"--> psurf = ",psurf
     405      ENDIF
     406      !! we need reference pressures
     407      pa=psurf/30.
     408      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
     409
     410!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     411!!!! END PLANETARY CONSTANTS !!!!
     412!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    258413
    259414c  Date et heure locale du debut du run
    … …  
    298453      ndt=ndt*day_step     
    299454      dtphys=daysec/day_step 
    300 
     455      write(*,*)"-------------------------------------"
     456      write(*,*)"-------------------------------------"
     457      write(*,*)"--> Physical timestep is ",dtphys
     458      write(*,*)"-------------------------------------"
     459      write(*,*)"-------------------------------------"
    301460
    302461c output spectrum?
    … …  
    306465      write(*,*)" specOLR = ",specOLR
    307466
    308 
    309 c Pression de surface sur la planete
    310 c ------------------------------------
    311 c
    312       psurf=100000. ! default value for psurf
    313       PRINT *,'Surface pressure (Pa) ?'
    314       call getin("psurf",psurf)
    315       write(*,*) " psurf = ",psurf
    316 c Pression de reference
    317       pa=psurf/30.
    318       preff=psurf ! these values are not needed in 1D  (are you sure JL12)
    319 
    320 c Gravity of planet
    321 c -----------------
    322       g=10.0 ! default value for g
    323       PRINT *,'Gravity ?'
    324       call getin("g",g)
    325       write(*,*) " g = ",g
    326 
    327       mugaz=0.
    328       cpp=0.
    329 
    330       call su_gases
    331       call calc_cpp_mugaz
    332      
    333 
     467!!! AS12: what shall we do with this?
    334468c Proprietes des poussiere aerosol
    335469c --------------------------------
    … …  
    339473      write(*,*) " tauvis = ",tauvis
    340474 
    341 c  latitude/longitude
    342 c  ------------------
    343       lati(1)=0 ! default value for lati(1)
    344       PRINT *,'latitude (in degrees) ?'
    345       call getin("latitude",lati(1))
    346       write(*,*) " latitude = ",lati(1)
    347       lati(1)=lati(1)*pi/180.E+0
    348       long(1)=0.E+0
    349       long(1)=long(1)*pi/180.E+0
    350 
    351 c  periastron/apoastron
    352 c  --------------------
    353       periastr=227.94
    354       apoastr=227.94 ! default = mean martian values
    355       PRINT *,'periastron (in AU) ?'
    356       call getin("periastr",periastr)
    357       write(*,*) "periastron = ",periastr
    358       periastr=periastr*149.598 ! AU to Mkm
    359       PRINT *,'apoastron (in AU) ?'
    360       call getin("apoastr",apoastr)
    361       write(*,*) "apoastron = ",apoastr
    362       apoastr=apoastr*149.598 ! AU to Mkm
    363 
    364 c  obliquity
    365 c  ---------
    366       obliquit=0.0! default=zero in 1d
    367       PRINT *,'obliquity (in AU) ?'
    368       call getin("obliquit",obliquit)
    369       write(*,*) "obliquity = ",obliquit
    370 
    371 c  Initialisation albedo / inertie du sol
    372 c  --------------------------------------
    373       albedodat(1)=0.2 ! default value for albedodat
    374       PRINT *,'Albedo of bare ground ?'
    375       call getin("albedo",albedodat(1))
    376       write(*,*) " albedo = ",albedodat(1)
    377 
    378       inertiedat(1,1)=400 ! default value for inertiedat
    379       PRINT *,'Soil thermal inertia (SI) ?'
    380       call getin("inertia",inertiedat(1,1))
    381       write(*,*) " inertia = ",inertiedat(1,1)
    382475c
    383476c  pour le schema d'ondes de gravite
    384477c  ---------------------------------
    385 c
    386478      zmea(1)=0.E+0
    387479      zstd(1)=0.E+0
    … …  
    403495      ENDDO
    404496
    405 c   Initialisation speciales "physiq"
    406 c   ---------------------------------
    407 c   la surface de chaque maille est inutile en 1D --->
    408       area(1)=1.E+0
    409       aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
    410 
    411 c   le geopotentiel au sol est inutile en 1D car tout est controle
    412 c   par la pression de surface --->
    413       phisfi(1)=0.E+0
    414 
    415 c  "inifis" reproduit un certain nombre d'initialisations deja faites
    416 c  + lecture des clefs de callphys.def
    417 c  + calcul de la frequence d'appel au rayonnement
    418 c  + calcul des sinus et cosinus des longitude latitude
    419 
    420 !Mars possible issue with dtphys in input and include!!!
    421       CALL inifis(1,llm,day0,daysec,dtphys,
    422      .            lati,long,area,rad,g,r,cpp)
    423497c   Initialisation pour prendre en compte les vents en 1-D
    424498c   ------------------------------------------------------
    … …  
    488562      write(*,*) " autozlevs = ", autozlevs
    489563
    490       pceil=100.0 ! Pascals
     564      pceil = psurf / 1000.0 ! Pascals
    491565      PRINT *,'Ceiling pressure (Pa) ?'
    492566      call getin("pceil",pceil)
    … …  
    571645     
    572646
     647c  Initialisation albedo / inertie du sol
     648c  --------------------------------------
     649      albedodat(1)=0.2 ! default value for albedodat
     650      PRINT *,'Albedo of bare ground ?'
     651      call getin("albedo",albedodat(1))
     652      write(*,*) " albedo = ",albedodat(1)
     653
     654      inertiedat(1,1)=400 ! default value for inertiedat
     655      PRINT *,'Soil thermal inertia (SI) ?'
     656      call getin("inertia",inertiedat(1,1))
     657      write(*,*) " inertia = ",inertiedat(1,1)
    573658
    574659! Initialize soil properties and temperature
Note: See TracChangeset for help on using the changeset viewer.