Changeset 3045 for trunk


Ignore:
Timestamp:
Sep 21, 2023, 10:19:40 AM (14 months ago)
Author:
jbclement
Message:

Mars PCM 1D:
Addition of the pressure profile to be able to run chained simulations in 1D. Like this, the surface pressure can be got from one run to the next.
Rework of "write_profile.F90" + correction of arguments in its call in "testphys1d.F": 'qsurf' was wrongly given!
JBC

Location:
trunk/LMDZ.MARS
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3042 r3045  
    41974197== 13/09/2023 == JBC
    41984198The date and time are redefined at the end of PEM to write the file "restarfi_evol.nc". Like so, the next GCM runs can restart at the date and time given at the beginning of the simulation since the PEM is only progressing by step of one Martian year.
     4199
     4200== 21/09/2023 == JBC
     4201Addition of the pressure profile to be able to run chained simulations in 1D. Like this, the surface pressure can be got from one run to the next.
     4202Rework of "write_profile.F90" + correction of arguments in its call in "testphys1d.F": 'qsurf' was wrongly given!

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F

    r3042 r3045  
    215215c  Prescribed constants to be set here
    216216c ------------------------------------------------------
    217 
    218217c      if(.not. startfile_1D ) then
    219 
    220218      pi=2.E+0*asin(1.E+0)
    221219
    … …  
    261259c     ----------------------------------------------------
    262260      cell_area(1)=1.E+0
    263 
    264261
    265262! check if there is a 'traceur.def' file
    … …  
    354351       WRITE(*,*) 'nqfils=',nqfils
    355352
    356         ! initialize tracers here:
     353        ! Initialize tracers here:
    357354        write(*,*) "testphys1d: initializing tracers"
    358         call read_profile(nq, nlayer, qsurf, q)
     355        call read_profile(nq,nlayer,qsurf,q)
    359356
    360357      call init_physics_distribution(regular_lonlat,4,
    … …  
    398395c  Discretization (Definition of grid and time steps)
    399396c  --------------
    400 c
    401397      nlevel=nlayer+1
    402398      nsoil=nsoilmx
    403399
    404400      day_step=48 ! default value for day_step
    405       PRINT *,'Number of time steps per sol ?'
     401      write(*,*)'Number of time steps per sol ?'
    406402      call getin("day_step",day_step)
    407403      write(*,*) " day_step = ",day_step
    … …  
    410406
    411407      ndt=10 ! default value for ndt
    412       PRINT *,'Number of sols to run ?'
     408      write(*,*)'Number of sols to run ?'
    413409      call getin("ndt",ndt)
    414410      write(*,*) " ndt = ",ndt
    … …  
    420416c Imposed surface pressure
    421417c ------------------------------------
    422 c
    423 
    424       psurf=610. ! default value for psurf
    425       PRINT *,'Surface pressure (Pa) ?'
    426       call getin("psurf",psurf)
    427       write(*,*) " psurf = ",psurf
     418      psurf = 610. ! Default value for psurf
     419      write(*,*) 'Surface pressure (Pa)?'
     420      open(3,file = 'profile_pressure',status = 'old',
     421     &     form = 'formatted',iostat = ierr)
     422      if (ierr == 0) then
     423          write(*,*) 'Reading file ''profile_pressure''...'
     424          read(3,*) psurf
     425          write(*,*) " psurf = ",psurf
     426      else
     427          write(*,*) 'File ''profile_pressure'' not found!'
     428          write(*,*) 'Reading surface pressure in ''callphys.def''',
     429     &' if given...'
     430          call getin("psurf",psurf)
     431          write(*,*) " psurf = ",psurf
     432      endif
     433      close(3)
    428434
    429435c Reference pressures
    … …  
    440446c Orbital parameters
    441447c ------------------
    442 
    443448      if(.not. startfile_1D ) then
    444449
    … …  
    449454        write(*,*) "Orbital parameters from callphys.def"
    450455        write(*,*) "Enter eccentricity & Lsperi"
    451         print *, 'Martian eccentricity (0<e<1) ?'
     456        write(*,*) 'Martian eccentricity (0<e<1) ?'
    452457        call getin('excentric ',excentric)
    453458        write(*,*)"excentric =",excentric
    454         print *, 'Solar longitude of perihelion (0<Ls<360) ?'
     459        write(*,*) 'Solar longitude of perihelion (0<Ls<360) ?'
    455460        call getin('Lsperi',Lsperi )
    456461        write(*,*)"Lsperi=",Lsperi
    … …  
    467472        write(*,*) "Default present-day orbital parameters"
    468473        write(*,*) "Unless specified otherwise"
    469         print *,'Min. distance Sun-Mars (Mkm)?'
     474        write(*,*)'Min. distance Sun-Mars (Mkm)?'
    470475        call getin("periheli",periheli)
    471476        write(*,*) " periheli = ",periheli
    472477
    473         print *,'Max. distance Sun-Mars (Mkm)?'
     478        write(*,*)'Max. distance Sun-Mars (Mkm)?'
    474479        call getin("aphelie",aphelie)
    475480        write(*,*) " aphelie = ",aphelie
    476481
    477         print *,'Day of perihelion?'
     482        write(*,*)'Day of perihelion?'
    478483        call getin("periday",peri_day)
    479484        write(*,*) " periday = ",peri_day
    480485
    481         print *,'Obliquity?'
     486        write(*,*)'Obliquity?'
    482487        call getin("obliquit",obliquit)
    483488        write(*,*) " obliquit = ",obliquit
    … …  
    489494c  ------------------
    490495      latitude(1)=0 ! default value for latitude
    491       PRINT *,'latitude (in degrees) ?'
     496      write(*,*)'latitude (in degrees) ?'
    492497      call getin("latitude",latitude(1))
    493498      write(*,*) " latitude = ",latitude
    … …  
    533538
    534539      albedodat(1)=0.2 ! default value for albedodat
    535       PRINT *,'Albedo of bare ground ?'
     540      write(*,*)'Albedo of bare ground ?'
    536541      call getin("albedo",albedodat(1))
    537542      write(*,*) " albedo = ",albedodat(1)
    … …  
    539544
    540545      inertiedat(1,1)=400 ! default value for inertiedat
    541       PRINT *,'Soil thermal inertia (SI) ?'
     546      write(*,*)'Soil thermal inertia (SI) ?'
    542547      call getin("inertia",inertiedat(1,1))
    543548      write(*,*) " inertia = ",inertiedat(1,1)
    … …  
    581586c 
    582587      hmons(1)=0.E+0
    583       PRINT *,'hmons is initialized to ',hmons(1)
     588      write(*,*)'hmons is initialized to ',hmons(1)
    584589      summit(1)=0.E+0
    585       PRINT *,'summit is initialized to ',summit(1)
     590      write(*,*)'summit is initialized to ',summit(1)
    586591      base(1)=0.E+0
    587592c
    … …  
    604609c    geostrophic wind
    605610      gru=10. ! default value for gru
    606       PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
     611      write(*,*)'zonal eastward component of the geostrophic',
     612     &' wind (m/s) ?'
    607613      call getin("u",gru)
    608614      write(*,*) " u = ",gru
    609615      grv=0. !default value for grv
    610       PRINT *,'meridional northward component of the geostrophic',
     616      write(*,*)'meridional northward component of the geostrophic',
    611617     &' wind (m/s) ?'
    612618      call getin("v",grv)
    … …  
    641647      if(.not. startfile_1D ) then
    642648      qsurf(igcm_co2)=0.E+0 ! default value for co2ice
    643       PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
     649      write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
    644650      call getin("co2ice",qsurf(igcm_co2))
    645651      write(*,*) " co2ice = ",qsurf(igcm_co2)
    … …  
    664670c    """"""""""""""""""""
    665671      hybrid=.true.
    666       PRINT *,'Hybrid coordinates ?'
     672      write(*,*)'Hybrid coordinates ?'
    667673      call getin("hybrid",hybrid)
    668674      write(*,*) " hybrid = ", hybrid
    … …  
    807813      endif !(.not. startfile_1D )
    808814      watercaptag(1)=.false. ! Default: no water ice reservoir
    809       print *,'Water ice cap on ground ?'
     815      write(*,*)'Water ice cap on ground ?'
    810816      call getin("watercaptag",watercaptag)
    811817      write(*,*) " watercaptag = ",watercaptag
    … …  
    816822      atm_wat_profile = -1. ! Default: free atm wat profile
    817823      if (water) then
    818       print *,'Force atmospheric water vapor profile?'
     824      write(*,*)'Force atmospheric water vapor profile?'
    819825      call getin("atm_wat_profile",atm_wat_profile)
    820826      write(*,*) "atm_wat_profile = ", atm_wat_profile
    … …  
    993999      ENDDO   ! of idt=1,ndt ! end of time stepping loop
    9941000
    995       ! update the profiles files at the end of the run
    996       write_prof=.false.
     1001      ! Update the profile files at the end of the run
     1002      write_prof = .false.
    9971003      call getin("write_prof",write_prof)
    998       IF (write_prof) then
    999               DO iq = 1,nq
    1000                 call writeprofile(nlayer,q(:,iq),noms(iq),qsurf)
    1001               ENDDO
    1002       ENDIF
     1004      if (write_prof) then
     1005          do iq = 1,nq
     1006              call writeprofile(nlayer,noms(iq),qsurf(iq),q(:,iq))
     1007          enddo
     1008          call writeprofile(nlayer,'pressure',psurf,play)
     1009      endif
    10031010
    10041011      write(*,*) "testphys1d: Everything is cool."

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/write_profile.F90

    r3000 r3045  
    1 SUBROUTINE writeprofile(nlev,data,name2,qsurf)
     1SUBROUTINE writeprofile(nlev,profilename,surfdata,profiledata)
    22
     3implicit none
    34
    4 IMPLICIT NONE
     5! arguments
     6integer, intent(in)               :: nlev
     7real, intent(in)                  :: surfdata
     8real, dimension(nlev), intent(in) :: profiledata
     9character(len = 30), intent(in)   :: profilename
    510
    6 ! arguments:
    7 INTEGER nlev
    8 REAL data(nlev)
    9 REAL qsurf
    10 CHARACTER(len=30) :: name2
     11! Local
     12integer :: il
    1113
    12 !local
    13 INTEGER il
     14! Write the data
     15open(1,file = 'profile_out_'//trim(profilename),form = 'formatted')
     16write(1,*) surfdata
     17do il = 1,nlev
     18    write(1,*) profiledata(il)
     19enddo
     20close(1)
     21return
    1422
    15 !write the data
    16 OPEN(1,file='profile_out_'//trim(name2),form='formatted')
    17 write(1,*) qsurf
    18 DO il=1,nlev
    19     write(1,*) data(il)
    20 ENDDO
    21 CLOSE(1)
    22 RETURN
    2323END
Note: See TracChangeset for help on using the changeset viewer.