Changeset 999 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Jul 2, 2013, 9:40:28 AM (12 years ago)

Author:

tnavarro

Message:

Possibility to store multiple initial states in one start/startfi. This is RETROCOMPATIBLE. New option ecrithist in run.def to write data in start/startfi every ecrithist dynamical timestep. New option timestart in run.def to initialize the GCM with the time timestart stored in start

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• iniwrite.F (modified) (1 diff)
• phyetat0.F (modified) (15 diffs)
• physdem.F (added)
• physiq.F (modified) (5 diffs)
• soil_settings.F (modified) (5 diffs)
• testphys1d.F (modified) (1 diff)

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

@@ -72 +72 @@
       tab_cntrl(15) = stot0
       tab_cntrl(16) = ang0
+
+      tab_cntrl(27) = hour_ini
 c
 c    ..........    P.Le Van  ( ajout le 8/04/96 )    .........

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

@@ -1 +1 @@
       SUBROUTINE phyetat0 (fichnom,tab0,Lmodif,nsoil,nq,
-     .           day_ini,time,
+     .           day_ini,time0,
      .           tsurf,tsoil,emis,q2,qsurf,co2ice)
 
@@ -13 +13 @@
 !                      r4 or r8 restarts independently of having compiled
 !                      the GCM in r4 or r8)
+!         June 2013 TN : Possibility to read files with a time axis
 c======================================================================
 !#include "netcdf.inc"
@@ -24 +25 @@
 #include "comcstfi.h"
 !#include "tracer.h"
-#include"advtrac.h"
+#include "advtrac.h"
+#include "control.h"
 c======================================================================
       INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
@@ -38 +40 @@
       integer nq
       integer day_ini
-      real time
+      real time0
 
 !  outputs:
@@ -72 +74 @@
       character(len=30) :: txt ! to store some text
 
- 
+! specific for time 
+      REAL,ALLOCATABLE :: time(:) ! times stored in start
+      INTEGER timelen ! number of times stored in the file
+      INTEGER indextime ! index of selected time
+
+      INTEGER edges(3),corner(3)
+
+
 c
 c Ouvrir le fichier contenant l etat initial:
@@ -111 +120 @@
       write(*,*) 'TABFI de phyeta0',Lmodif,tab0
       call tabfi (nid,Lmodif,tab0,day_ini,lmax,p_rad,
-     .              p_omeg,p_g,p_mugaz,p_daysec,time)
+     .              p_omeg,p_g,p_mugaz,p_daysec,time0)
 c
 c Read latitudes (coordinates): No need, these are provided by the dynamics
@@ -333 +342 @@
       ENDDO
       PRINT*,'<zthe>:', xmin, xmax
+      
+c
+c Time axis
+c
+      ierr = nf90_inq_dimid(nid,"Time",nvarid)
+
+      IF (ierr .NE. nf90_noerr) THEN
+        indextime = 1
+        print*, "phyetat0: No time axis found in "//fichnom
+      
+      ELSE
+        print*, "phyetat0: Time axis found in "//fichnom
+       
+        ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
+        
+        ALLOCATE(time(timelen))
+  
+        ierr=nf90_inq_varid(nid,"Time",nvarid)
+        IF (ierr .NE. nf90_noerr) THEN
+             ierr=nf90_inq_varid(nid,"temps", nvarid)
+                 IF (ierr .NE. nf90_noerr) THEN
+             PRINT*, "dynetat0: <Time> or <temps> absent"
+             write(*,*)trim(nf90_strerror(ierr))
+             CALL abort
+           ENDIF
+        ENDIF
+        ierr=nf90_get_var(nid,nvarid,time)
+        IF (ierr .NE. nf90_noerr) THEN
+           PRINT*, "dynetat0: Lecture echouee <Time>/<temps>"
+           write(*,*)trim(nf90_strerror(ierr))
+           CALL abort
+        ENDIF
+c      
+c Select desired time
+c
+        IF (timestart .lt. 0) THEN  ! default: we use the last time value
+          indextime = timelen
+        ELSE  ! else we look for the desired value in the time axis
+         indextime = 0
+          DO i=1,timelen
+            IF (abs(time(i) - timestart) .lt. 0.01) THEN
+               indextime = i
+               EXIT
+            ENDIF
+          ENDDO
+          IF (indextime .eq. 0) THEN
+            PRINT*, "Time", timestart," is not in "//fichnom//"!!"
+            PRINT*, "Stored times are:"
+            DO i=1,timelen
+               PRINT*, time(i)
+            ENDDO
+            CALL abort
+          ENDIF
+        ENDIF
+        
+        ! In startfi the absolute date is day_ini + time0 + time
+        ! For now on, in the GCM physics, it is day_ini + time0
+        time0 = time(indextime) + time0
+        day_ini = day_ini + INT(time0)
+        time0 = time0 - INT(time0)
+        
+        PRINT*, "phyetat0: Selected time ",time(indextime),
+     .          " at index ",indextime
+     
+        DEALLOCATE(time)
+        
+      ENDIF ! of IF (Time not in file)
+      
+      
 c
 c CO2 ice cover
 c
+      corner(1)=1
+      corner(2)=indextime
+      edges(1)=ngridmx
+      edges(2)=1
       ierr=nf90_inq_varid(nid,"co2ice",nvarid)
       IF (ierr.NE.nf90_noerr) THEN
@@ -342 +424 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,co2ice)
+      ierr=nf90_get_var(nid,nvarid,co2ice,corner,edges)
       IF (ierr.NE.nf90_noerr) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <co2ice>'
@@ -389 +471 @@
          PRINT*, '          J ignore donc les autres temperatures TS**'
          !ierr=nf90_get_var(nid,nvarid,tsurf(1,1))
-         ierr=nf90_get_var(nid,nvarid,tsurf)
+         ierr=nf90_get_var(nid,nvarid,tsurf,corner,edges)
          IF (ierr.NE.nf90_noerr) THEN
             PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
@@ -453 +535 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,emis)
+      ierr=nf90_get_var(nid,nvarid,emis,corner,edges)
       IF (ierr.NE.nf90_noerr) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <emis>'
@@ -491 +573 @@
 c pbl wind variance
 c
+      corner(1)=1
+      corner(2)=1
+      corner(3)=indextime
+      edges(1)=ngridmx
+      edges(2)=llm+1
+      edges(3)=1
       ierr=nf90_inq_varid(nid,"q2",nvarid)
       IF (ierr.NE.nf90_noerr) THEN
@@ -497 +585 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,q2)
+      ierr=nf90_get_var(nid,nvarid,q2,corner,edges)
       IF (ierr.NE.nf90_noerr) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <q2>'
@@ -511 +599 @@
 c tracer on surface
 c
-
+      corner(1)=1
+      corner(2)=indextime
+      edges(1)=ngridmx
+      edges(2)=1
       IF(nq.GE.1) THEN
          nqold=nq
@@ -542 +633 @@
            ELSE
            !ierr=nf90_get_var(nid,nvarid,qsurf(1,iq))
-           ierr=nf90_get_var(nid,nvarid,surffield)
+           ierr=nf90_get_var(nid,nvarid,surffield,corner,edges)
            qsurf(1:ngridmx,iq)=surffield(1:ngridmx)
              IF (ierr.NE.nf90_noerr) THEN
@@ -583 +674 @@
 ! as well as thermal inertia and volumetric heat capacity
 
-      call soil_settings(nid,ngridmx,nsoil,tsurf,tsoil)
+      call soil_settings(nid,ngridmx,nsoil,tsurf,tsoil,indextime)
 c
 c Fermer le fichier:

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

@@ -162 +162 @@
 
       REAL pday
-      REAL ptime 
+      REAL ptime
       logical tracerdyn
 
@@ -214 +214 @@
       REAL sky
 
-      SAVE day_ini, icount
+      SAVE day_ini, icount, time_phys
       SAVE aerosol, tsurf,tsoil
       SAVE co2ice,albedo,emis, q2
@@ -478 +478 @@
 #ifndef MESOSCALE
          if (callslope) call getslopes(phisfi)
+                           
+         call physdem0("restartfi.nc",long,lati,nsoilmx,nq,
+     .                 ptimestep,pday,time_phys,area,
+     .                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+     
 #endif
                    
@@ -1354 +1359 @@
       
 
+
 c-----------------------------------------------------------------------
 c  10. Write output files
@@ -1520 +1526 @@
 c              thus we store for time=time+dtvr
 
-         IF(lastcall) THEN
-            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 
-            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
-     .              ptimestep,pday,
-     .              ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
-     .              area,albedodat,inertiedat,zmea,zstd,zsig,
-     .              zgam,zthe)
+         IF(   ((ecritstart.GT.0) .and. 
+     .          (MOD(icount*iphysiq,ecritstart).EQ.0)) 
+     .    .or. lastcall  ) THEN
+         
+          ztime_fin = pday + ptime  + ptimestep/(float(iphysiq)*daysec)
+     .               - day_ini - time_phys
+          print*, pday,ptime,day_ini, time_phys
+          write(*,'(A,I7,A,F12.5)') 
+     .         'PHYSIQ: Ecriture du fichier restartfi ; icount=',
+     .          icount,' date=',ztime_fin
+            
+            
+          call physdem1("restartfi.nc",nsoilmx,nq,
+     .                ptimestep,ztime_fin,
+     .                tsurf,tsoil,co2ice,emis,q2,qsurf)
+          
          ENDIF
 #endif

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

@@ -1 @@
-      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
+      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
 
       use netcdf
@@ -12 +12 @@
 !                      r4 or r8 restarts independently of having compiled
 !                      the GCM in r4 or r8)
+!                June 2013 TN : Possibility to read files with a time axis
+!
 !
 !  This subroutine reads from a NetCDF file (opened by the caller)
@@ -42 +44 @@
       integer nsoil     ! # of soil layers
       real tsurf(ngrid) ! surface temperature
+      integer indextime ! position on time axis
 !  output:
       real tsoil(ngridmx,nsoilmx)       ! soil temperature
@@ -53 +56 @@
       integer ndims     ! # of dimensions of read <inertiedat> data
       integer ig,iloop  ! loop counters
+      
+      integer edges(3),corner(3) ! to read a specific time
 
       logical :: olddepthdef=.false. ! flag
@@ -298 +303 @@
         endif
        else ! put values in tsoil
+        corner(1)=1
+        corner(2)=1
+        corner(3)=indextime
+        edges(1)=ngridmx
+        edges(2)=nsoilmx
+        edges(3)=1
+        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
         ierr=nf90_get_var(nid,nvarid,tsoil)
         if (ierr.ne.nf90_noerr) then

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

@@ -641 +641 @@
 c  It is needed to transfert physics variables to "physiq"...
 
-      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
-     .              dtphys,float(day0),time,tsurf,
-     .              tsoil,co2ice,emis,q2,qsurf,area,albedodat,
-     .              inertiedat,zmea,zstd,zsig,zgam,zthe)
+      call physdem0("startfi.nc",long,lati,nsoilmx,nqmx,
+     .              dtphys,float(day0),time,area,
+     .              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+      call physdem1("startfi.nc",nsoilmx,nqmx,
+     .              dtphys,time,
+     .              tsurf,tsoil,co2ice,emis,q2,qsurf)
 
 c=======================================================================