Changeset 3051

Timestamp:

Sep 26, 2023, 6:10:40 PM (16 months ago)

Author:

jbclement

Message:

Mars PCM 1D:
Addition of the file "start1D.txt" which mimics a "start.nc" file to be able to run chained simulations in 1D. Like this, key variables and profiles can be got from one run to the next. As a consequence, the subroutine "write_profile.F90" is replaced by "writerestart1D.F90".
JBC

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/dyn1d/testphys1d.F (modified) (33 diffs)
• libf/phymars/dyn1d/writerestart1D.F90 (moved) (moved from trunk/LMDZ.MARS/libf/phymars/dyn1d/write_profile.F90) (1 diff)
• libf/phymars/phyredem.F90 (modified) (1 diff)

trunk/LMDZ.MARS/changelog.txt

@@ -4207 +4207 @@
 
 == 26/09/2023 == JBC
+Mars PEM:
 Minor changes concerning the form of the code in the PEM.
+Mars PCM 1D:
+Addition of the file "start1D.txt" which mimics a "start.nc" file to be able to run chained simulations in 1D. Like this, key variables and profiles can be got from one run to the next. As a consequence, the subroutine "write_profile.F90" is replaced by "writerestart1D.F90".

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

@@ -90 +90 @@
       INTEGER ilayer,ilevel,isoil,idt,iq
       LOGICAl firstcall,lastcall
-      LOGICAL write_prof
 c
       real,parameter :: odpref=610. ! DOD reference pressure (Pa)
@@ -140 +139 @@
       REAL halfaxe, excentric, Lsperi
       Logical paleomars
-c   JN : Force atmospheric water profiles
-      REAL atm_wat_profile
-      REAL atm_wat_tau
-      REAL,ALLOCATABLE :: zqsat(:) ! useful for (atm_wat_profile=2)
-
+
+c   JN & JBC: Force atmospheric water profiles
+      real                            :: atm_wat_profile, atm_wat_tau
+      real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
 
 c   MVals: isotopes as in the dynamics (CRisi)
@@ -152 +150 @@
       INTEGER ierr0
       LOGICAL :: continu
-
       logical :: present
 
@@ -158 +155 @@
       real :: flux_geo_tmp
 
-c   RV: Start from a startfi and not run.def
-      logical :: startfile_1D
-      REAL tab_cntrl(100)
-      LOGICAL :: found
-      REAL albedo_read(1,2,1)      ! surface albedo
+c   RV & JBC: Use of "start1D.txt" and "startfi.nc" files
+      logical                :: startfiles_1D, found
+      logical                :: therestart1D, therestartfi
+      character(len = 30)    :: header
+      real, dimension(100)   :: tab_cntrl
+      real, dimension(1,2,1) :: albedo_read(1,2,1) ! surface albedo
 
 c   LL: Possibility to add subsurface ice
@@ -184 +182 @@
 !      call initcomgeomphy
 
-
 c ------------------------------------------------------
 c  Loading run parameters from "run.def" file
 c ------------------------------------------------------
-
-
 ! check if 'run.def' file is around (otherwise reading parameters
 ! from callphys.def via getin() routine won't work.
@@ -206 +201 @@
       endif
 
-      write(*,*)'Do you want to start with a startfi and write
-     & a restartfi?'
-      startfile_1D=.false.
-      call getin("startfile_1D",startfile_1D)
-      write(*,*)" startfile_1D = ", startfile_1D
+      write(*,*)'Do you want to use "start1D.txt" and "startfi.nc"',
+     &' files?'
+      startfiles_1D = .false.
+      call getin("startfiles_1D",startfiles_1D)
+      write(*,*) " startfiles_1D = ", startfiles_1D
+      
+      if (startfiles_1D) then
+          inquire(file ='start1D.txt',exist = therestart1D) 
+          if (.not. therestart1D) then
+              write(*,*) 'There is no "start1D.txt" file!'
+              write(*,*) 'Initialization is done with default values.'
+          endif
+          inquire(file ='startfi.nc',exist = therestartfi) 
+          if (.not. therestartfi) then
+              write(*,*) 'There is no "startfi.nc" file!'
+              write(*,*) 'Initialization is done with default values.'
+          endif
+      endif
 
 c ------------------------------------------------------
 c  Prescribed constants to be set here
 c ------------------------------------------------------
-c      if(.not. startfile_1D ) then
       pi=2.E+0*asin(1.E+0)
 
@@ -253 +260 @@
       dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
       dtemisice(2) = 2.          ! time scale for snow metamorphism (south)
-
-computeTcondc      endif ! .not. startfile_1D
 
 c     mesh surface (not a very usefull quantity in 1D)
@@ -353 +358 @@
         ! Initialize tracers here:
         write(*,*) "testphys1d: initializing tracers"
-        call read_profile(nq,nlayer,qsurf,q)
+        if (.not. therestart1D) then
+            call read_profile(nq,nlayer,qsurf,q)
+        else
+            do iq = 1,nq
+                open(3,file = 'start1D.txt',status = "old",
+     &action = "read")
+                read(3,*) header, qsurf(iq),
+     & (q(ilayer,iq), ilayer = 1,nlayer)
+                if (tname(iq) /= trim(header)) then
+                    write(*,*) 'Tracer names not compatible for',
+     &' initialization with "start1D.txt"!'
+                    stop
+                endif
+            enddo
+        endif
 
       call init_physics_distribution(regular_lonlat,4,
      &                               1,1,1,nlayer,COMM_LMDZ)
 
-      if(.not. startfile_1D ) then
-
 c  Date and local time at beginning of run
 c  ---------------------------------------
+      if (.not. startfiles_1D) then
 c    Date (in sols since spring solstice) at beginning of run
       day0 = 0 ! default value for day0
@@ -391 +409 @@
        call close_startphy
 
-      endif !startfile_1D
+      endif !startfiles_1D
 
 c  Discretization (Definition of grid and time steps)
@@ -418 +436 @@
       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
+      if (.not. therestart1D) then
+          call getin("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
+          read(3,*) header, psurf
       endif
-      close(3)
+      write(*,*) " psurf = ",psurf
 
 c Reference pressures
@@ -446 +456 @@
 c Orbital parameters
 c ------------------
-      if(.not. startfile_1D ) then
-
+      if (.not. startfiles_1D) then
       paleomars=.false. ! Default: no water ice reservoir
       call getin("paleomars",paleomars)
@@ -487 +496 @@
         call getin("obliquit",obliquit)
         write(*,*) " obliquit = ",obliquit
-      end if
-
-      endif !(.not. startfile_1D )
+      endif
+
+      endif !(.not. startfiles_1D )
  
 c  latitude/longitude
@@ -533 +542 @@
 c  Initialize albedo / soil thermal inertia
 c  ----------------------------------------
-c
-
-      if(.not. startfile_1D ) then
-
+      if (.not. startfiles_1D) then
       albedodat(1)=0.2 ! default value for albedodat
       write(*,*)'Albedo of bare ground ?'
@@ -556 +562 @@
       write(*,*) " z0 = ",z0(1)
 
-      endif !(.not. startfile_1D ) 
+      endif !(.not. startfiles_1D ) 
 
 ! Initialize local slope parameters (only matters if "callslope"
@@ -575 +581 @@
 c  for the gravity wave scheme
 c  ---------------------------------
-c
       zmea(1)=0.E+0
       zstd(1)=0.E+0
@@ -584 +589 @@
 c  for the slope wind scheme
 c  ---------------------------------
-c  
       hmons(1)=0.E+0
       write(*,*)'hmons is initialized to ',hmons(1)
@@ -593 +597 @@
 c  Default values initializing the coefficients calculated later
 c  ---------------------------------
-c  
       tauscaling(1)=1. ! calculated in aeropacity_mod.F
       totcloudfrac(1)=1. ! calculated in watercloud_mod.F      
@@ -619 +622 @@
       write(*,*) " v = ",grv
 
-c     Initialize winds  for first time step
-      DO ilayer=1,nlayer
-         u(ilayer)=gru
-         v(ilayer)=grv
-         w(ilayer)=0 ! default: no vertical wind
-      ENDDO
+c     Initialize winds for first time step
+      if (.not. therestart1D) then
+          do ilayer = 1,nlayer
+             u(ilayer) = gru
+             v(ilayer) = grv
+          enddo
+      else
+          read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
+          read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
+      endif
+      w = 0. ! default: no vertical wind
 
 c     Initialize turbulente kinetic energy
@@ -645 +653 @@
       endif
 
-      if(.not. startfile_1D ) then
+      if (.not. startfiles_1D) then
       qsurf(igcm_co2)=0.E+0 ! default value for co2ice
       write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
       call getin("co2ice",qsurf(igcm_co2))
       write(*,*) " co2ice = ",qsurf(igcm_co2)
-      endif !(.not. startfile_1D )
+      endif !(.not. startfiles_1D )
 
 c
 c  emissivity
 c  ----------
-      if(.not. startfile_1D ) then
+      if (.not. startfiles_1D) then
       emis=emissiv
       IF (qsurf(igcm_co2).eq.1.E+0) THEN
@@ -661 +669 @@
          IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
       ENDIF
-      endif !(.not. startfile_1D ) 
+      endif !(.not. startfiles_1D ) 
  
-
 c  Compute pressures and altitudes of atmospheric levels 
 c  ----------------------------------------------------------------
-
 c    Vertical Coordinates
 c    """"""""""""""""""""
@@ -705 +711 @@
       call profile(nlayer+1,tmp1,tmp2)
 
-      tsurf=tmp2(0)
-      DO ilayer=1,nlayer
-        temp(ilayer)=tmp2(ilayer)
-      ENDDO
-      
+      if (.not. therestart1D) then
+          tsurf = tmp2(0)
+          do ilayer = 1,nlayer
+              temp(ilayer) = tmp2(ilayer)
+          enddo
+      else
+          read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer)
+          close(3)
+      endif      
 
 ! Initialize soil properties and temperature
@@ -715 +725 @@
       volcapa=1.e6 ! volumetric heat capacity
 
-      if(.not. startfile_1D ) then
-
+      if (.not. startfiles_1D) then
 ! Initialize depths
 ! -----------------
@@ -766 +775 @@
        ENDDO
 
-      endif !(.not. startfile_1D )
+      endif !(.not. startfiles_1D)
 
       flux_geo_tmp=0.
@@ -781 +790 @@
       enddo
 
-c    Initialize traceurs
-c    ---------------------------
-
+c Initialize traceurs
+c ---------------------------
       if (photochem.or.callthermos) then
          write(*,*) 'Initializing chemical species'
@@ -809 +817 @@
 c Check if the surface is a water ice reservoir 
 c --------------------------------------------------
-      if(.not. startfile_1D ) then
-      watercap(1,:)=0 ! Initialize watercap
-      endif !(.not. startfile_1D )
+      if (.not. startfiles_1D) watercap(1,:)=0 ! Initialize watercap
       watercaptag(1)=.false. ! Default: no water ice reservoir
       write(*,*)'Water ice cap on ground ?'
@@ -822 +828 @@
       atm_wat_profile = -1. ! Default: free atm wat profile
       if (water) then
-      write(*,*)'Force atmospheric water vapor profile?'
-      call getin("atm_wat_profile",atm_wat_profile)
-      write(*,*) "atm_wat_profile = ", atm_wat_profile
-      if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
-        write(*,*) "Free atmospheric water vapor profile"
-        write(*,*) "Total water is conserved in the column"
-      else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
-        write(*,*) "Dry atmospheric water vapor profile"
-      else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
-        write(*,*) "Prescribed atmospheric water vapor profile"
-        write(*,*) "Unless it reaches saturation (maximal value)"
-      else
-        write(*,*) 'Water vapor profile value not correct!'
-        stop
-      endif
+          write(*,*)'Force atmospheric water vapor profile?'
+          call getin('atm_wat_profile',atm_wat_profile)
+          write(*,*) 'atm_wat_profile = ', atm_wat_profile
+          if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
+              write(*,*) 'Free atmospheric water vapor profile'
+              write(*,*) 'Total water is conserved in the column'
+          else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
+              write(*,*) 'Dry atmospheric water vapor profile'
+          else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.)
+     & then
+              write(*,*) 'Prescribed atmospheric water vapor profile'
+              write(*,*) 'Unless it reaches saturation (maximal value)'
+          else
+              write(*,*) 'Water vapor profile value not correct!'
+              stop
+          endif
       endif
 
@@ -844 +851 @@
       atm_wat_tau = -1. ! Default: no time relaxation
       if (water) then
-      write(*,*) 'Relax atmospheric water vapor profile?'
-      call getin("atm_wat_tau",atm_wat_tau)
-      write(*,*) "atm_wat_tau = ", atm_wat_tau
-      if (atm_wat_tau < 0.) then
-        write(*,*) "Atmospheric water vapor profile is not relaxed"
-      else
-        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
-          write(*,*) "Relaxed atmospheric water vapor profile towards ",
-     &    atm_wat_profile
-          write(*,*) "Unless it reaches saturation (maximal value)"
-        else
-       write(*,*) 'Reference atmospheric water vapor profile not known!'
-       write(*,*) 'Please, specify atm_wat_profile'
-       stop
-        endif
-      endif
+          write(*,*) 'Relax atmospheric water vapor profile?'
+          call getin('atm_wat_tau',atm_wat_tau)
+          write(*,*) 'atm_wat_tau = ', atm_wat_tau
+          if (atm_wat_tau < 0.) then
+              write(*,*) 'Atmospheric water vapor profile is not',
+     &' relaxed.'
+          else
+              if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.)
+     & then
+                  write(*,*) 'Relaxed atmospheric water vapor profile',
+     &' towards ', atm_wat_profile
+                  write(*,*) 'Unless it reaches saturation (maximal',
+     &' value)'
+              else
+                  write(*,*) 'Reference atmospheric water vapor',
+     &' profile not known!'
+                  write(*,*) 'Please, specify atm_wat_profile'
+                  stop
+              endif
+          endif
       endif
 
@@ -867 +878 @@
 c  It is needed to transfert physics variables to "physiq"...
 
-      if(.not. startfile_1D ) then
-
+      if (.not. startfiles_1D) then
       call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,
      &              llm,nq,dttestphys,float(day0),0.,cell_area,
@@ -877 +887 @@
      &              q2,qsurf,tauscaling,
      &              totcloudfrac,wstar,watercap,perenial_co2ice)
-      endif !(.not. startfile_1D )
+      endif !(.not. startfiles_1D )
 
 c=======================================================================
@@ -919 +929 @@
 !       ---------------------------------------
       if (water) then
-        call watersat(nlayer,temp,play,zqsat)
-        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
-        ! If atmospheric water is monitored
-          if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile
-            ! Wet if >0, Dry if =0
-            q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf)
-            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
-          else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
-        q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap)
-     &          - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau)
-            q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap))
-            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
+          call watersat(nlayer,temp,play,zqsat)
+          if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
+          ! If atmospheric water is monitored
+              if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile: wet if >0, dry if =0
+              q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf)
+              q(:,igcm_h2o_ice) = 0. ! reset h2o ice
+              else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
+       q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap)
+     &        - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau)
+              q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap))
+              q(:,igcm_h2o_ice) = 0. ! reset h2o ice
+              endif
           endif
-        endif
       endif
 
@@ -939 +948 @@
 c       --------------------
       CALL physiq (1,llm,nq,firstcall,lastcall,
-     ,     day,time,dttestphys,plev,play,phi,
-     ,     u,v,temp,q,w,
+     .     day,time,dttestphys,plev,play,phi,
+     .     u,v,temp,q,w,
 C - outputs
-     s     du, dv, dtemp, dq,dpsurf)
+     .     du, dv, dtemp, dq,dpsurf)
 !      write(*,*) "testphys1d apres q", q(1,:)
 
@@ -992 +1001 @@
 
 !       increment tracers
-        DO iq = 1, nq
+        DO iq = 1,nq
           DO ilayer=1,nlayer
              q(ilayer,iq)=q(ilayer,iq)+dttestphys*dq(ilayer,iq)
@@ -999 +1008 @@
       ENDDO   ! of idt=1,ndt ! end of time stepping loop
 
-      ! 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,noms(iq),qsurf(iq),q(:,iq))
-          enddo
-          call writeprofile(nlayer,'pressure',psurf,play)
-      endif
+      ! Writing the "restart1D.txt" file for the next run
+      if (startfiles_1D) call writerestart1D(psurf,tsurf,nlayer,temp,
+     &u,v,nq,noms,qsurf,q)
+ 
 
       write(*,*) "testphys1d: Everything is cool."

trunk/LMDZ.MARS/libf/phymars/dyn1d/writerestart1D.F90

@@ -1 @@
-SUBROUTINE writeprofile(nlev,profilename,surfdata,profiledata)
+SUBROUTINE writerestart1D(psurf,tsurf,nlayer,temp,u,v,nq,qnames,qsurf,q)
 
 implicit none
 
-! arguments
-integer, intent(in)               :: nlev
-real, intent(in)                  :: surfdata
-real, dimension(nlev), intent(in) :: profiledata
-character(len = 30), intent(in)   :: profilename
+! Arguments
+integer, intent(in)                           :: nlayer, nq
+real, intent(in)                              :: psurf, tsurf
+real, dimension(nlayer), intent(in)           :: temp, u, v
+real, dimension(nlayer,nq), intent(in)        :: q
+real, dimension(nq), intent(in)               :: qsurf
+character(len = *), dimension(nq), intent(in) :: qnames
 
-! Local
-integer :: il
+! Local variables
+integer :: il, iq
 
-! Write the data
-open(1,file = 'profile_out_'//trim(profilename),form = 'formatted')
-write(1,*) surfdata
-do il = 1,nlev
-    write(1,*) profiledata(il)
+! Write the data needed for a restart in "restart1D.txt"
+open(1,file = 'restart1D.txt',status = "replace",action = "write")
+do iq = 1,nq
+    write(1,*) qnames(iq), qsurf(iq), (q(il,iq), il = 1,nlayer)
 enddo
+write(1,*) 'ps', psurf
+write(1,*) 'u', (u(il), il = 1,nlayer)
+write(1,*) 'v', (v(il), il = 1,nlayer)
+write(1,*) 'teta', tsurf, (temp(il), il = 1,nlayer)
 close(1)
-return
 
-END
+END SUBROUTINE writerestart1D

trunk/LMDZ.MARS/libf/phymars/phyredem.F90

@@ -224 +224 @@
  
  ! Perenial CO2 ice layer
-  if(paleoclimate) then
-    call put_field("perenial_co2ice","CO2 ice cover",perenial_co2ice,time)
-  endif
+  if (paleoclimate) call put_field("perenial_co2ice","CO2 ice cover",perenial_co2ice,time)
+
   ! Surface temperature
   call put_field("tsurf","Surface temperature",tsurf,time)