Changeset 999 for trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F

Timestamp:

Jul 2, 2013, 9:40:28 AM (11 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

File:

• trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F (modified) (18 diffs)

trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F

@@ -1 +1 @@
       SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
-     .                    teta,q,masse,ps,phis,time)
+     .                    teta,q,masse,ps,phis,time0)
       
       use netcdf
@@ -20 +20 @@
 c                    r4 or r8 restarts independently of having compiled
 c                    the GCM in r4 or r8)
-c
+c           June 2013 TN : Possibility to read files with a time axis.
+c                          NOW defrun_new MUST BE CALLED BEFORE dynetat0.
 c=======================================================================
 c-----------------------------------------------------------------------
@@ -36 +37 @@
 #include "serre.h"
 #include "logic.h"
-#include"advtrac.h"
+#include "advtrac.h"
+#include "control.h"
 
 c   Arguments:
@@ -48 +50 @@
       REAL ps(iip1,jjp1),phis(iip1,jjp1)
 
-      REAL time ! fraction of day the fields correspond to
+      REAL time0 ! fraction of day the fields correspond to
 
 c   Variables 
@@ -58 +60 @@
       INTEGER ierr, nid, nvarid, nqold
       CHARACTER  str3*3,yes*1
+      
+      REAL,ALLOCATABLE :: time(:) ! times stored in start
+      INTEGER timelen ! number of times stored in the file
+      INTEGER indextime ! index of selected time
+      !REAL  hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1
+
+      INTEGER edges(4),corner(4)
 
 c-----------------------------------------------------------------------
@@ -102 +111 @@
       pa         = tab_cntrl(17)
       preff      = tab_cntrl(18)
+      
+      hour_ini   = tab_cntrl(29)
+
 c
       clon       = tab_cntrl(19)
@@ -119 +131 @@
         IF( tab_cntrl(26).EQ.1. ) ysinus = . TRUE. 
       ENDIF
+      
 c   .................................................................
 c
@@ -239 +252 @@
       ENDIF
 
+
+! read time axis
+      ierr = nf90_inq_dimid(nid,"Time",nvarid)
+      ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
+
+      ALLOCATE(time(timelen))
+
       ierr=nf90_inq_varid(nid,"Time",nvarid)
       IF (ierr .NE. nf90_noerr) THEN
@@ -254 +274 @@
          CALL abort
       ENDIF
-
+      
+! select desired time
+      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 start the absolute date is day_ini + hour_ini + time
+      ! For now on, in the GCMdynamics, it is day_ini + time0
+      time0 = time(indextime) + hour_ini
+      day_ini = day_ini + INT(time0)
+      time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
+      hour_ini = time0
+      
+      PRINT*, "dynetat0: Selected time ",time(indextime),
+     .        " at index ",indextime
+      
+      DEALLOCATE(time)
+
+
+! read vcov
+      corner(1)=1
+      corner(2)=1
+      corner(3)=1
+      corner(4)=indextime
+      edges(1)=iip1
+      edges(2)=jjm
+      edges(3)=llm
+      edges(4)=1
+      ierr=nf90_inq_varid(nid,"vcov",nvarid)
+      IF (ierr .NE. nf90_noerr) THEN
+         PRINT*, "dynetat0: Le champ <vcov> est absent"
+         write(*,*)trim(nf90_strerror(ierr))
+         CALL abort
+      ENDIF
+      !ierr=nf90_get_var(nid,nvarid,vcov)
+      ierr=nf90_get_var(nid,nvarid,vcov,corner,edges)
+      IF (ierr .NE. nf90_noerr) THEN
+         PRINT*, "dynetat0: Lecture echouee pour <vcov>"
+         write(*,*)trim(nf90_strerror(ierr))
+         CALL abort
+      ENDIF
+      
+! read ucov
+      corner(1)=1
+      corner(2)=1
+      corner(3)=1
+      corner(4)=indextime
+      edges(1)=iip1
+      edges(2)=jjp1
+      edges(3)=llm
+      edges(4)=1
       ierr=nf90_inq_varid(nid,"ucov",nvarid)
       IF (ierr .NE. nf90_noerr) THEN
@@ -261 +347 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,ucov)
+      ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
       IF (ierr .NE. nf90_noerr) THEN
          PRINT*, "dynetat0: Lecture echouee pour <ucov>"
@@ -267 +353 @@
          CALL abort
       ENDIF
- 
-      ierr=nf90_inq_varid(nid,"vcov",nvarid)
-      IF (ierr .NE. nf90_noerr) THEN
-         PRINT*, "dynetat0: Le champ <vcov> est absent"
-         write(*,*)trim(nf90_strerror(ierr))
-         CALL abort
-      ENDIF
-      ierr=nf90_get_var(nid,nvarid,vcov)
-      IF (ierr .NE. nf90_noerr) THEN
-         PRINT*, "dynetat0: Lecture echouee pour <vcov>"
-         write(*,*)trim(nf90_strerror(ierr))
-         CALL abort
-      ENDIF
-
+
+! read teta
       ierr=nf90_inq_varid(nid,"teta",nvarid)
       IF (ierr .NE. nf90_noerr) THEN
@@ -287 +361 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,teta)
+      ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
       IF (ierr .NE. nf90_noerr) THEN
          PRINT*, "dynetat0: Lecture echouee pour <teta>"
@@ -294 +368 @@
       ENDIF
 
-
+! read tracers
       IF(nq.GE.1) THEN
         write(*,*) 'dynetat0: loading tracers'
@@ -318 +392 @@
            ELSE
            !ierr=nf90_get_var(nid,nvarid,q(1,1,1,iq))
-           ierr=nf90_get_var(nid,nvarid,traceur)
+           ierr=nf90_get_var(nid,nvarid,traceur,corner,edges)
              IF (ierr .NE. nf90_noerr) THEN
 !                 PRINT*, "dynetat0: Lecture echouee pour "//str3
@@ -355 +429 @@
       ENDIF
 
+! read masse
       ierr=nf90_inq_varid(nid,"masse",nvarid)
       IF (ierr .NE. nf90_noerr) THEN
@@ -361 +436 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,masse)
+      ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
       IF (ierr .NE. nf90_noerr) THEN
          PRINT*, "dynetat0: Lecture echouee pour <masse>"
@@ -368 +443 @@
       ENDIF
 
+! read ps
+      corner(1)=1
+      corner(2)=1
+      corner(3)=indextime
+      edges(1)=iip1
+      edges(2)=jjp1
+      edges(3)=1
       ierr=nf90_inq_varid(nid,"ps",nvarid)
       IF (ierr .NE. nf90_noerr) THEN
@@ -374 +456 @@
          CALL abort
       ENDIF
-      ierr=nf90_get_var(nid,nvarid,ps)
+      ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
       IF (ierr .NE. nf90_noerr) THEN
          PRINT*, "dynetat0: Lecture echouee pour <ps>"