Changeset 3045 for trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F

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

File:

• trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F (modified) (19 diffs)

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

@@ -215 +215 @@
 c  Prescribed constants to be set here
 c ------------------------------------------------------
-
 c      if(.not. startfile_1D ) then
-
       pi=2.E+0*asin(1.E+0)
 
@@ -261 +259 @@
 c     ----------------------------------------------------
       cell_area(1)=1.E+0
-
 
 ! check if there is a 'traceur.def' file
@@ -354 +351 @@
        WRITE(*,*) 'nqfils=',nqfils
 
-        ! initialize tracers here:
+        ! Initialize tracers here:
         write(*,*) "testphys1d: initializing tracers"
-        call read_profile(nq, nlayer, qsurf, q)
+        call read_profile(nq,nlayer,qsurf,q)
 
       call init_physics_distribution(regular_lonlat,4,
@@ -398 +395 @@
 c  Discretization (Definition of grid and time steps)
 c  --------------
-c
       nlevel=nlayer+1
       nsoil=nsoilmx
 
       day_step=48 ! default value for day_step
-      PRINT *,'Number of time steps per sol ?'
+      write(*,*)'Number of time steps per sol ?'
       call getin("day_step",day_step)
       write(*,*) " day_step = ",day_step
@@ -410 +406 @@
 
       ndt=10 ! default value for ndt
-      PRINT *,'Number of sols to run ?'
+      write(*,*)'Number of sols to run ?'
       call getin("ndt",ndt)
       write(*,*) " ndt = ",ndt
@@ -420 +416 @@
 c Imposed surface pressure
 c ------------------------------------
-c
-
-      psurf=610. ! default value for psurf
-      PRINT *,'Surface pressure (Pa) ?'
-      call getin("psurf",psurf)
-      write(*,*) " psurf = ",psurf
+      psurf = 610. ! Default value for psurf
+      write(*,*) 'Surface pressure (Pa)?'
+      open(3,file = 'profile_pressure',status = 'old',
+     &     form = 'formatted',iostat = ierr)
+      if (ierr == 0) then
+          write(*,*) 'Reading file ''profile_pressure''...'
+          read(3,*) psurf
+          write(*,*) " psurf = ",psurf
+      else
+          write(*,*) 'File ''profile_pressure'' not found!'
+          write(*,*) 'Reading surface pressure in ''callphys.def''',
+     &' if given...'
+          call getin("psurf",psurf)
+          write(*,*) " psurf = ",psurf
+      endif
+      close(3)
 
 c Reference pressures
@@ -440 +446 @@
 c Orbital parameters
 c ------------------
-
       if(.not. startfile_1D ) then
 
@@ -449 +454 @@
         write(*,*) "Orbital parameters from callphys.def"
         write(*,*) "Enter eccentricity & Lsperi"
-        print *, 'Martian eccentricity (0<e<1) ?'
+        write(*,*) 'Martian eccentricity (0<e<1) ?'
         call getin('excentric ',excentric)
         write(*,*)"excentric =",excentric 
-        print *, 'Solar longitude of perihelion (0<Ls<360) ?'
+        write(*,*) 'Solar longitude of perihelion (0<Ls<360) ?'
         call getin('Lsperi',Lsperi )
         write(*,*)"Lsperi=",Lsperi 
@@ -467 +472 @@
         write(*,*) "Default present-day orbital parameters"
         write(*,*) "Unless specified otherwise"
-        print *,'Min. distance Sun-Mars (Mkm)?'
+        write(*,*)'Min. distance Sun-Mars (Mkm)?'
         call getin("periheli",periheli)
         write(*,*) " periheli = ",periheli
 
-        print *,'Max. distance Sun-Mars (Mkm)?'
+        write(*,*)'Max. distance Sun-Mars (Mkm)?'
         call getin("aphelie",aphelie)
         write(*,*) " aphelie = ",aphelie
 
-        print *,'Day of perihelion?'
+        write(*,*)'Day of perihelion?'
         call getin("periday",peri_day)
         write(*,*) " periday = ",peri_day
 
-        print *,'Obliquity?'
+        write(*,*)'Obliquity?'
         call getin("obliquit",obliquit)
         write(*,*) " obliquit = ",obliquit
@@ -489 +494 @@
 c  ------------------
       latitude(1)=0 ! default value for latitude
-      PRINT *,'latitude (in degrees) ?'
+      write(*,*)'latitude (in degrees) ?'
       call getin("latitude",latitude(1))
       write(*,*) " latitude = ",latitude
@@ -533 +538 @@
 
       albedodat(1)=0.2 ! default value for albedodat
-      PRINT *,'Albedo of bare ground ?'
+      write(*,*)'Albedo of bare ground ?'
       call getin("albedo",albedodat(1))
       write(*,*) " albedo = ",albedodat(1)
@@ -539 +544 @@
 
       inertiedat(1,1)=400 ! default value for inertiedat
-      PRINT *,'Soil thermal inertia (SI) ?'
+      write(*,*)'Soil thermal inertia (SI) ?'
       call getin("inertia",inertiedat(1,1))
       write(*,*) " inertia = ",inertiedat(1,1)
@@ -581 +586 @@
 c  
       hmons(1)=0.E+0
-      PRINT *,'hmons is initialized to ',hmons(1)
+      write(*,*)'hmons is initialized to ',hmons(1)
       summit(1)=0.E+0
-      PRINT *,'summit is initialized to ',summit(1)
+      write(*,*)'summit is initialized to ',summit(1)
       base(1)=0.E+0
 c
@@ -604 +609 @@
 c    geostrophic wind
       gru=10. ! default value for gru
-      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
+      write(*,*)'zonal eastward component of the geostrophic',
+     &' wind (m/s) ?'
       call getin("u",gru)
       write(*,*) " u = ",gru
       grv=0. !default value for grv
-      PRINT *,'meridional northward component of the geostrophic',
+      write(*,*)'meridional northward component of the geostrophic',
      &' wind (m/s) ?'
       call getin("v",grv)
@@ -641 +647 @@
       if(.not. startfile_1D ) then
       qsurf(igcm_co2)=0.E+0 ! default value for co2ice
-      PRINT *,'Initial CO2 ice on the surface (kg.m-2)'
+      write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
       call getin("co2ice",qsurf(igcm_co2))
       write(*,*) " co2ice = ",qsurf(igcm_co2)
@@ -664 +670 @@
 c    """"""""""""""""""""
       hybrid=.true.
-      PRINT *,'Hybrid coordinates ?'
+      write(*,*)'Hybrid coordinates ?'
       call getin("hybrid",hybrid)
       write(*,*) " hybrid = ", hybrid
@@ -807 +813 @@
       endif !(.not. startfile_1D )
       watercaptag(1)=.false. ! Default: no water ice reservoir
-      print *,'Water ice cap on ground ?'
+      write(*,*)'Water ice cap on ground ?'
       call getin("watercaptag",watercaptag)
       write(*,*) " watercaptag = ",watercaptag
@@ -816 +822 @@
       atm_wat_profile = -1. ! Default: free atm wat profile
       if (water) then
-      print *,'Force atmospheric water vapor profile?'
+      write(*,*)'Force atmospheric water vapor profile?'
       call getin("atm_wat_profile",atm_wat_profile)
       write(*,*) "atm_wat_profile = ", atm_wat_profile
@@ -993 +999 @@
       ENDDO   ! of idt=1,ndt ! end of time stepping loop
 
-      ! update the profiles files at the end of the run
-      write_prof=.false.
+      ! Update the profile files at the end of the run
+      write_prof = .false.
       call getin("write_prof",write_prof)
-      IF (write_prof) then
-              DO iq = 1,nq
-                call writeprofile(nlayer,q(:,iq),noms(iq),qsurf)
-              ENDDO
-      ENDIF
+      if (write_prof) then
+          do iq = 1,nq
+              call writeprofile(nlayer,noms(iq),qsurf(iq),q(:,iq))
+          enddo
+          call writeprofile(nlayer,'pressure',psurf,play)
+      endif
 
       write(*,*) "testphys1d: Everything is cool."