Changeset 2924 for trunk/LMDZ.MARS

Timestamp:

Mar 30, 2023, 10:37:08 AM (20 months ago)

Author:

romain.vande

Message:

1D Mars PCM:
Add the possibility to start the 1d model from a startfi (option startfile_1D in run.def false by default).
It will then write a restartfi at the end of the run.
+ Correction of r2919 for flux geo in 1D and initialize it to 0 by default.
+ Minor correction of iostart error message.
RV

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/dyn1d/testphys1d.F (modified) (22 diffs)
• libf/phymars/iostart.F90 (modified) (3 diffs)
• libf/phymars/physiq_mod.F (modified) (5 diffs)

trunk/LMDZ.MARS/README

@@ -3964 +3964 @@
 Add illustrative execution job for spirit (aka MESOIPSL cluster) in deftank
 
+==30/03/2023 == RV
+Add the possibility to start the 1d model from a startfi (option startfile_1D in run.def false by default).
+It will then write a restartfi at the end of the run.
++ Correction of r2919 for flux geo in 1D and initialize it to 0 by default.
++ Minor correction of iostart error message.
+

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

@@ -6 +6 @@
       use mod_grid_phy_lmdz, only : regular_lonlat
       use infotrac, only: nqtot, tname, nqperes,nqfils
-      use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx,flux_geo
+      use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx,
+     &                     flux_geo
       use comgeomfi_h, only: sinlat, ini_fillgeom
       use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
@@ -13 +14 @@
      &                     watercaptag, watercap, hmons, summit, base
       use slope_mod, only: theta_sl, psi_sl
-      use comslope_mod, only: def_slope,subslope_dist
+      use comslope_mod, only: def_slope,subslope_dist,def_slope_mean
       use phyredem, only: physdem0,physdem1
       use geometry_mod, only: init_geometry
@@ -36 +37 @@
       USE read_profile_mod, ONLY: read_profile
       use inichim_newstart_mod, only: inichim_newstart
+      use phyetat0_mod, only: phyetat0
+      USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE,
+     &                  nf90_strerror
+      use iostart, only: open_startphy,get_var, close_startphy
       IMPLICIT NONE
 
@@ -146 +151 @@
 c   LL: Subsurface geothermal flux
       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=======================================================================
 
@@ -156 +167 @@
 !      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.
+      open(99,file='run.def',status='old',form='formatted',
+     &     iostat=ierr)
+      if (ierr.ne.0) then
+        write(*,*) 'Cannot find required file "run.def"'
+        write(*,*) '  (which should contain some input parameters'
+        write(*,*) '   along with the following line:'
+        write(*,*) '   INCLUDEDEF=callphys.def'
+        write(*,*) '   )'
+        write(*,*) ' ... might as well stop here ...'
+        stop
+      else
+        close(99)
+      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
+
 c ------------------------------------------------------
 c  Prescribed constants to be set here
 c ------------------------------------------------------
+
+c      if(.not. startfile_1D ) then
 
       pi=2.E+0*asin(1.E+0)
@@ -196 +237 @@
       iceradius(2) = 100.e-6     ! mean scat radius of CO2 snow (south)
       dtemisice(1) = 2.          ! time scale for snow metamorphism (north)
-      dtemisice(2) = 2.          ! time scale for snow metamorphism (south
+      dtemisice(2) = 2.          ! time scale for snow metamorphism (south)
+
+c      endif ! .not. startfile_1D
 
 c     mesh surface (not a very usefull quantity in 1D)
 c     ----------------------------------------------------
       cell_area(1)=1.E+0
-      
-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.
-      open(99,file='run.def',status='old',form='formatted',
-     &     iostat=ierr)
-      if (ierr.ne.0) then
-        write(*,*) 'Cannot find required file "run.def"'
-        write(*,*) '  (which should contain some input parameters'
-        write(*,*) '   along with the following line:'
-        write(*,*) '   INCLUDEDEF=callphys.def'
-        write(*,*) '   )'
-        write(*,*) ' ... might as well stop here ...'
-        stop
-      else
-        close(99)
-      endif
 
 
@@ -320 +342 @@
         call read_profile(nq, nlayer, qsurf, q)
 
+      call init_physics_distribution(regular_lonlat,4,
+     &                               1,1,1,nlayer,1)
+
+      if(.not. startfile_1D ) then
+
 c  Date and local time at beginning of run
 c  ---------------------------------------
@@ -335 +362 @@
       time=time/24.E+0 ! convert time (hours) to fraction of sol
 
+      else
+      call open_startphy("startfi.nc")
+      call get_var("controle",tab_cntrl,found)
+       if (.not.found) then
+         call abort_physic("open_startphy", 
+     &        "tabfi: Failed reading <controle> array",1)
+       else
+         write(*,*)'tabfi: tab_cntrl',tab_cntrl
+       endif
+       day0 = tab_cntrl(3)
+       day=float(day0)
+
+       call get_var("Time",time,found)
+
+       call close_startphy
+
+      endif !startfile_1D
+
 c  Discretization (Definition of grid and time steps)
 c  --------------
@@ -360 +405 @@
 c ------------------------------------
 c
+
       psurf=610. ! default value for psurf
       PRINT *,'Surface pressure (Pa) ?'
       call getin("psurf",psurf)
       write(*,*) " psurf = ",psurf
+
 c Reference pressures
       pa=20.   ! transition pressure (for hybrid coord.)
@@ -377 +424 @@
 c Orbital parameters
 c ------------------
+
+      if(.not. startfile_1D ) then
+
       paleomars=.false. ! Default: no water ice reservoir
       call getin("paleomars",paleomars)
@@ -417 +467 @@
         write(*,*) " obliquit = ",obliquit
       end if
+
+      endif !(.not. startfile_1D )
  
 c  latitude/longitude
@@ -433 +485 @@
 !Mars possible matter with dtphys in input and include!!!
 ! Initializations below should mimick what is done in iniphysiq for 3D GCM
-      call init_physics_distribution(regular_lonlat,4,
-     &                               1,1,1,nlayer,1)
       call init_interface_dyn_phys
       call init_regular_lonlat(1,1,longitude,latitude,
@@ -463 +513 @@
 c  ----------------------------------------
 c
+
+      if(.not. startfile_1D ) then
+
       albedodat(1)=0.2 ! default value for albedodat
       PRINT *,'Albedo of bare ground ?'
@@ -478 +531 @@
       call getin("z0",z0(1))
       write(*,*) " z0 = ",z0(1)
+
+      endif !(.not. startfile_1D ) 
 
 ! Initialize local slope parameters (only matters if "callslope"
@@ -564 +619 @@
         stop
       endif
+
+      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)'
       call getin("co2ice",qsurf(igcm_co2))
       write(*,*) " co2ice = ",qsurf(igcm_co2)
+      endif !(.not. startfile_1D )
 
 c
 c  emissivity
 c  ----------
+      if(.not. startfile_1D ) then
       emis=emissiv
       IF (qsurf(igcm_co2).eq.1.E+0) THEN
@@ -577 +636 @@
          IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere
       ENDIF
-
+      endif !(.not. startfile_1D ) 
  
 
@@ -632 +691 @@
 ! ------------------------------------------
       volcapa=1.e6 ! volumetric heat capacity
+      if(.not. startfile_1D ) then
       DO isoil=1,nsoil
          inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
          tsoil(isoil)=tsurf(1)  ! soil temperature
       ENDDO
-
+      endif !(.not. startfile_1D )
+
+      flux_geo_tmp=0.
       call getin("flux_geo",flux_geo_tmp)
       flux_geo(:,:) = flux_geo_tmp
 
-
-
-
-flux_geo
 ! Initialize depths
 ! -----------------
@@ -681 +739 @@
 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 )
       watercaptag(1)=.false. ! Default: no water ice reservoir
       print *,'Water ice cap on ground ?'
@@ -723 +783 @@
 c  It is needed to transfert physics variables to "physiq"...
 
-      call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
-     &              nq,dtphys,float(day0),0.,cell_area,
+      if(.not. startfile_1D ) then
+
+      call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,
+     &              llm,nq,dtphys,float(day0),0.,cell_area,
      &              albedodat,inertiedat,def_slope,subslope_dist)
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
@@ -730 +792 @@
      &              tsurf,tsoil,albedo,emis,q2,qsurf,tauscaling,
      &              totcloudfrac,wstar,watercap)
+      else
+         CALL phyetat0 ("startfi.nc",0,0,
+     &         nsoilmx,ngrid,nlayer,nq,
+     &         day0,time,
+     &         tsurf,tsoil,albedo_read,emis,
+     &         q2,qsurf,tauscaling,totcloudfrac,wstar,
+     &         watercap,def_slope,def_slope_mean,subslope_dist)
+         day=real(day0)
+         albedo(1,1)=albedo_read(1,1,1)
+      endif !(.not. startfile_1D )
 
 c=======================================================================
@@ -737 +809 @@
       firstcall=.true.
       lastcall=.false.
-
       DO idt=1,ndt
-c        IF (idt.eq.ndt) lastcall=.true.
-        IF (idt.eq.ndt-day_step-1) then       !test
-         lastcall=.true.
-         call solarlong(day*1.0,zls)
-         write(103,*) 'Ls=',zls*180./pi
-         write(103,*) 'Lat=', latitude(1)*180./pi
-         write(103,*) 'Tau=', tauvis/odpref*psurf
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(103,*) 'RunEnd - Atmos. Temp. File'
-         write(104,*) 'Ls=',zls*180./pi
-         write(104,*) 'Lat=', latitude(1)
-         write(104,*) 'Tau=', tauvis/odpref*psurf
-         write(104,*) 'RunEnd - Atmos. Temp. File'
-        ENDIF
+        IF (idt.eq.ndt) lastcall=.true.
+c        IF (idt.eq.ndt-day_step-1) then       !test
+c         lastcall=.true.
+c         call solarlong(day*1.0,zls)
+c         write(103,*) 'Ls=',zls*180./pi
+c         write(103,*) 'Lat=', latitude(1)*180./pi
+c         write(103,*) 'Tau=', tauvis/odpref*psurf
+c         write(103,*) 'RunEnd - Atmos. Temp. File'
+c         write(103,*) 'RunEnd - Atmos. Temp. File'
+c         write(104,*) 'Ls=',zls*180./pi
+c         write(104,*) 'Lat=', latitude(1)
+c         write(104,*) 'Tau=', tauvis/odpref*psurf
+c         write(104,*) 'RunEnd - Atmos. Temp. File'
+c        ENDIF
 
 c     compute geopotential
@@ -879 +950 @@
         CALL endg1d(1,nlayer,zlay/1000.,ndt)
 
+c      if(startfile_1D) then
+c       call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid,
+c     &              llm,nq,dtphys,float(day0),0.,cell_area,
+c     &              albedodat,inertiedat,def_slope,subslope_dist)
+c       call physdem1("restartfi.nc",nsoilmx,ngrid,llm,nq,
+c     &              dtphys,time,
+c     &              tsurf,tsoil,albedo,emis,q2,qsurf,tauscaling,
+c     &              totcloudfrac,wstar,watercap)
+c      endif! startfile_1D
+
       write(*,*) "testphys1d: Everything is cool."
 

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

@@ -424 +424 @@
   
       ierr=NF90_INQ_VARID(nid_start,var_name,varid)
-      
       IF (ierr==NF90_NOERR) THEN
         ierr=NF90_GET_VAR(nid_start,varid,var)
@@ -1049 +1048 @@
       ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,(/idim1d/),nvarid)
 #endif
+      
       ! Add a "title" attribute
       IF (LEN_TRIM(title)>0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
@@ -1058 +1058 @@
         write(*,*)'put_var_rgen: problem writing '//trim(var_name)
         write(*,*)trim(nf90_strerror(ierr))
-        CALL abort_physic("get_var_rgen","Failed writing variable",1)
+        CALL abort_physic("put_var_rgen","Failed writing variable",1)
       ENDIF
     ENDIF ! of IF (is_master)

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -109 +109 @@
      &                        major_slope,compute_meshgridavg,
      &                        ini_comslope_h
-
+      USE ioipsl_getincom, only: getin
       IMPLICIT NONE
 c=======================================================================
@@ -545 +545 @@
 
       logical :: write_restart
+      logical,save :: startfile_1D
 
 c=======================================================================
@@ -748 +749 @@
 
 #ifndef MESOSCALE
-
-         if (ngrid.ne.1) then
+      startfile_1D=.false.
+      call getin("startfile_1D",startfile_1D)
+         if (ngrid.ne.1 .or. startfile_1D) then
            ! no need to compute slopes when in 1D; it is an input
            if (callslope) call getslopes(ngrid,phisfi)
@@ -2549 +2551 @@
 c        ----------------------------------------------------------
 
-      IF (ngrid.NE.1) THEN
+        IF (startfile_1D) THEN
 
 #ifndef MESOSCALE
@@ -2610 +2612 @@
      .                emis,q2,qsurf,tauscaling,totcloudfrac,wstar,
      .                watercap)
-          
-         ENDIF ! of IF (write_restart)
+          ENDIF ! of IF (write_restart)
+
 #endif
+
+         ENDIF ! startfile_1D
+
+         IF (ngrid.NE.1) then
 
 c        -------------------------------------------------------------------