source: trunk/LMDZ.MARS/libf/dyn3d/dynetat0.F @ 999

Last change on this file since 999 was 999, checked in by tnavarro, 11 years ago

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 size: 15.0 KB
Line 
1      SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
2     .                    teta,q,masse,ps,phis,time0)
3     
4      use netcdf
5     
6      IMPLICIT NONE
7
8c=======================================================================
9c
10c   Auteur:  P. Le Van / L.Fairhead
11c   -------
12c
13c   objet:
14c   ------
15c
16c   Lecture de l'etat initial
17c
18c   Modifs: Oct.2008 read in tracers by name. Ehouarn Millour
19c           Aug.2010 use NetCDF90 to load variables (enables using
20c                    r4 or r8 restarts independently of having compiled
21c                    the GCM in r4 or r8)
22c           June 2013 TN : Possibility to read files with a time axis.
23c                          NOW defrun_new MUST BE CALLED BEFORE dynetat0.
24c=======================================================================
25c-----------------------------------------------------------------------
26c   Declarations:
27c   -------------
28
29#include "dimensions.h"
30#include "paramet.h"
31#include "temps.h"
32#include "comconst.h"
33#include "comvert.h"
34#include "comgeom2.h"
35#include "ener.h"
36#include "description.h"
37#include "serre.h"
38#include "logic.h"
39#include "advtrac.h"
40#include "control.h"
41
42c   Arguments:
43c   ----------
44
45      CHARACTER*(*) fichnom
46      INTEGER nq
47      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm),teta(iip1,jjp1,llm)
48      REAL q(iip1,jjp1,llm,nq)
49      REAL masse(iip1,jjp1,llm)
50      REAL ps(iip1,jjp1),phis(iip1,jjp1)
51
52      REAL time0 ! fraction of day the fields correspond to
53
54c   Variables
55c
56      REAL traceur(iip1,jjp1,llm) ! to temporarily store a tracer
57      INTEGER length,iq,i,j,l
58      PARAMETER (length = 100)
59      REAL tab_cntrl(length) ! tableau des parametres du run
60      INTEGER ierr, nid, nvarid, nqold
61      CHARACTER  str3*3,yes*1
62     
63      REAL,ALLOCATABLE :: time(:) ! times stored in start
64      INTEGER timelen ! number of times stored in the file
65      INTEGER indextime ! index of selected time
66      !REAL  hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1
67
68      INTEGER edges(4),corner(4)
69
70c-----------------------------------------------------------------------
71
72c  Ouverture NetCDF du fichier etat initial
73
74      ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
75      IF (ierr.NE.nf90_noerr) THEN
76        write(6,*)' Pb d''ouverture du fichier ',trim(fichnom)
77        write(*,*)trim(nf90_strerror(ierr))
78        CALL ABORT
79      ENDIF
80
81c
82      ierr=nf90_inq_varid(nid,"controle",nvarid)
83      IF (ierr .NE. nf90_noerr) THEN
84         PRINT*, "dynetat0: Le champ <controle> est absent"
85         write(*,*)trim(nf90_strerror(ierr))
86         CALL abort
87      ENDIF
88      ierr=nf90_get_var(nid,nvarid,tab_cntrl)
89      IF (ierr .NE. nf90_noerr) THEN
90         PRINT*, "dynetat0: Lecture echoue pour <controle>"
91         write(*,*)trim(nf90_strerror(ierr))
92         CALL abort
93      ENDIF
94
95      im         = tab_cntrl(1)
96      jm         = tab_cntrl(2)
97      lllm       = tab_cntrl(3)
98      day_ini    = tab_cntrl(4)
99      rad        = tab_cntrl(5)
100      omeg       = tab_cntrl(6)
101      g          = tab_cntrl(7)
102      cpp        = tab_cntrl(8)
103      kappa      = tab_cntrl(9)
104      daysec     = tab_cntrl(10)
105      dtvr       = tab_cntrl(11)
106      etot0      = tab_cntrl(12)
107      ptot0      = tab_cntrl(13)
108      ztot0      = tab_cntrl(14)
109      stot0      = tab_cntrl(15)
110      ang0       = tab_cntrl(16)
111      pa         = tab_cntrl(17)
112      preff      = tab_cntrl(18)
113     
114      hour_ini   = tab_cntrl(29)
115
116c
117      clon       = tab_cntrl(19)
118      clat       = tab_cntrl(20)
119      grossismx  = tab_cntrl(21)
120      grossismy  = tab_cntrl(22)
121c
122      IF ( tab_cntrl(23).EQ.1. )  THEN
123        fxyhypb  = . TRUE .
124        dzoomx   = tab_cntrl(24)
125        dzoomy   = tab_cntrl(25)
126        taux     = tab_cntrl(27)
127        tauy     = tab_cntrl(28)
128      ELSE
129        fxyhypb = . FALSE .
130        ysinus  = . FALSE .
131        IF( tab_cntrl(26).EQ.1. ) ysinus = . TRUE.
132      ENDIF
133     
134c   .................................................................
135c
136c
137      PRINT*,'dynetat0: rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
138 
139      IF(   im.ne.iim           )  THEN
140          PRINT 1,im,iim
141          STOP
142      ELSE  IF( jm.ne.jjm       )  THEN
143          PRINT 2,jm,jjm
144          STOP
145      ELSE  IF( lllm.ne.llm     )  THEN
146          PRINT 3,lllm,llm
147          STOP
148      ENDIF
149
150      ierr=nf90_inq_varid(nid,"rlonu",nvarid)
151      IF (ierr .NE. nf90_noerr) THEN
152         PRINT*, "dynetat0: Le champ <rlonu> est absent"
153         write(*,*)trim(nf90_strerror(ierr))
154         CALL abort
155      ENDIF
156      ierr=nf90_get_var(nid,nvarid,rlonu)
157      IF (ierr .NE. nf90_noerr) THEN
158         PRINT*, "dynetat0: Lecture echouee pour <rlonu>"
159         write(*,*)trim(nf90_strerror(ierr))
160         CALL abort
161      ENDIF
162
163      ierr=nf90_inq_varid(nid,"rlatu",nvarid)
164      IF (ierr .NE. nf90_noerr) THEN
165         PRINT*, "dynetat0: Le champ <rlatu> est absent"
166         write(*,*)trim(nf90_strerror(ierr))
167         CALL abort
168      ENDIF
169      ierr=nf90_get_var(nid,nvarid,rlatu)
170      IF (ierr .NE. nf90_noerr) THEN
171         PRINT*, "dynetat0: Lecture echouee pour <rlatu>"
172         write(*,*)trim(nf90_strerror(ierr))
173         CALL abort
174      ENDIF
175
176      ierr=nf90_inq_varid(nid,"rlonv",nvarid)
177      IF (ierr .NE. nf90_noerr) THEN
178         PRINT*, "dynetat0: Le champ <rlonv> est absent"
179         write(*,*)trim(nf90_strerror(ierr))
180         CALL abort
181      ENDIF
182      ierr=nf90_get_var(nid,nvarid,rlonv)
183      IF (ierr .NE. nf90_noerr) THEN
184         PRINT*, "dynetat0: Lecture echouee pour <rlonv>"
185         write(*,*)trim(nf90_strerror(ierr))
186         CALL abort
187      ENDIF
188
189      ierr=nf90_inq_varid(nid,"rlatv",nvarid)
190      IF (ierr .NE. nf90_noerr) THEN
191         PRINT*, "dynetat0: Le champ <rlatv> est absent"
192         write(*,*)trim(nf90_strerror(ierr))
193         CALL abort
194      ENDIF
195      ierr=nf90_get_var(nid,nvarid,rlatv)
196      IF (ierr .NE. nf90_noerr) THEN
197         PRINT*, "dynetat0: Lecture echouee pour rlatv"
198         write(*,*)trim(nf90_strerror(ierr))
199         CALL abort
200      ENDIF
201
202      ierr=nf90_inq_varid(nid,"cu",nvarid)
203      IF (ierr .NE. nf90_noerr) THEN
204         PRINT*, "dynetat0: Le champ <cu> est absent"
205         write(*,*)trim(nf90_strerror(ierr))
206         CALL abort
207      ENDIF
208      ierr=nf90_get_var(nid,nvarid,cu)
209      IF (ierr .NE. nf90_noerr) THEN
210         PRINT*, "dynetat0: Lecture echouee pour <cu>"
211         write(*,*)trim(nf90_strerror(ierr))
212         CALL abort
213      ENDIF
214
215      ierr=nf90_inq_varid(nid,"cv",nvarid)
216      IF (ierr .NE. nf90_noerr) THEN
217         PRINT*, "dynetat0: Le champ <cv> est absent"
218         write(*,*)trim(nf90_strerror(ierr))
219         CALL abort
220      ENDIF
221      ierr=nf90_get_var(nid,nvarid,cv)
222      IF (ierr .NE. nf90_noerr) THEN
223         PRINT*, "dynetat0: Lecture echouee pour <cv>"
224         write(*,*)trim(nf90_strerror(ierr))
225         CALL abort
226      ENDIF
227
228      ierr=nf90_inq_varid(nid,"aire",nvarid)
229      IF (ierr .NE. nf90_noerr) THEN
230         PRINT*, "dynetat0: Le champ <aire> est absent"
231         write(*,*)trim(nf90_strerror(ierr))
232         CALL abort
233      ENDIF
234      ierr=nf90_get_var(nid,nvarid,aire)
235      IF (ierr .NE. nf90_noerr) THEN
236         PRINT*, "dynetat0: Lecture echouee pour <aire>"
237         write(*,*)trim(nf90_strerror(ierr))
238         CALL abort
239      ENDIF
240
241      ierr=nf90_inq_varid(nid,"phisinit",nvarid)
242      IF (ierr .NE. nf90_noerr) THEN
243         PRINT*, "dynetat0: Le champ <phisinit> est absent"
244         write(*,*)trim(nf90_strerror(ierr))
245         CALL abort
246      ENDIF
247      ierr=nf90_get_var(nid,nvarid,phis)
248      IF (ierr .NE. nf90_noerr) THEN
249         PRINT*, "dynetat0: Lecture echouee pour <phisinit>"
250         write(*,*)trim(nf90_strerror(ierr))
251         CALL abort
252      ENDIF
253
254
255! read time axis
256      ierr = nf90_inq_dimid(nid,"Time",nvarid)
257      ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
258
259      ALLOCATE(time(timelen))
260
261      ierr=nf90_inq_varid(nid,"Time",nvarid)
262      IF (ierr .NE. nf90_noerr) THEN
263             ierr=nf90_inq_varid(nid,"temps", nvarid)
264                 IF (ierr .NE. nf90_noerr) THEN
265           PRINT*, "dynetat0: <Time> or <temps> absent"
266           write(*,*)trim(nf90_strerror(ierr))
267           CALL abort
268         ENDIF
269      ENDIF
270      ierr=nf90_get_var(nid,nvarid,time)
271      IF (ierr .NE. nf90_noerr) THEN
272         PRINT*, "dynetat0: Lecture echouee <Time>/<temps>"
273         write(*,*)trim(nf90_strerror(ierr))
274         CALL abort
275      ENDIF
276     
277! select desired time
278      IF (timestart .lt. 0) THEN  ! default: we use the last time value
279        indextime = timelen
280      ELSE  ! else we look for the desired value in the time axis
281       indextime = 0
282        DO i=1,timelen
283          IF (abs(time(i) - timestart) .lt. 0.01) THEN
284             indextime = i
285             EXIT
286          ENDIF
287        ENDDO
288        IF (indextime .eq. 0) THEN
289          PRINT*, "Time", timestart," is not in "//fichnom//"!!"
290          PRINT*, "Stored times are:"
291          DO i=1,timelen
292             PRINT*, time(i)
293          ENDDO
294          CALL abort
295        ENDIF
296      ENDIF
297     
298      ! In start the absolute date is day_ini + hour_ini + time
299      ! For now on, in the GCMdynamics, it is day_ini + time0
300      time0 = time(indextime) + hour_ini
301      day_ini = day_ini + INT(time0)
302      time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
303      hour_ini = time0
304     
305      PRINT*, "dynetat0: Selected time ",time(indextime),
306     .        " at index ",indextime
307     
308      DEALLOCATE(time)
309
310
311! read vcov
312      corner(1)=1
313      corner(2)=1
314      corner(3)=1
315      corner(4)=indextime
316      edges(1)=iip1
317      edges(2)=jjm
318      edges(3)=llm
319      edges(4)=1
320      ierr=nf90_inq_varid(nid,"vcov",nvarid)
321      IF (ierr .NE. nf90_noerr) THEN
322         PRINT*, "dynetat0: Le champ <vcov> est absent"
323         write(*,*)trim(nf90_strerror(ierr))
324         CALL abort
325      ENDIF
326      !ierr=nf90_get_var(nid,nvarid,vcov)
327      ierr=nf90_get_var(nid,nvarid,vcov,corner,edges)
328      IF (ierr .NE. nf90_noerr) THEN
329         PRINT*, "dynetat0: Lecture echouee pour <vcov>"
330         write(*,*)trim(nf90_strerror(ierr))
331         CALL abort
332      ENDIF
333     
334! read ucov
335      corner(1)=1
336      corner(2)=1
337      corner(3)=1
338      corner(4)=indextime
339      edges(1)=iip1
340      edges(2)=jjp1
341      edges(3)=llm
342      edges(4)=1
343      ierr=nf90_inq_varid(nid,"ucov",nvarid)
344      IF (ierr .NE. nf90_noerr) THEN
345         PRINT*, "dynetat0: Le champ <ucov> est absent"
346         write(*,*)trim(nf90_strerror(ierr))
347         CALL abort
348      ENDIF
349      ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
350      IF (ierr .NE. nf90_noerr) THEN
351         PRINT*, "dynetat0: Lecture echouee pour <ucov>"
352         write(*,*)trim(nf90_strerror(ierr))
353         CALL abort
354      ENDIF
355
356! read teta
357      ierr=nf90_inq_varid(nid,"teta",nvarid)
358      IF (ierr .NE. nf90_noerr) THEN
359         PRINT*, "dynetat0: Le champ <teta> est absent"
360         write(*,*)trim(nf90_strerror(ierr))
361         CALL abort
362      ENDIF
363      ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
364      IF (ierr .NE. nf90_noerr) THEN
365         PRINT*, "dynetat0: Lecture echouee pour <teta>"
366         write(*,*)trim(nf90_strerror(ierr))
367         CALL abort
368      ENDIF
369
370! read tracers
371      IF(nq.GE.1) THEN
372        write(*,*) 'dynetat0: loading tracers'
373         IF(nq.GT.99) THEN
374            PRINT*, "Trop de traceurs"
375            CALL abort
376         ENDIF
377         nqold=nq
378         DO iq=1,nq
379!           str3(1:1)='q'
380!           WRITE(str3(2:3),'(i2.2)') iq
381!           ierr =  NF_INQ_VARID (nid, str3, nvarid)
382! NB: tracers are now read in using their name ('tnom' from advtrac.h)
383!           write(*,*) "  loading tracer:",trim(tnom(iq))
384           ierr=nf90_inq_varid(nid,tnom(iq),nvarid)
385           IF (ierr .NE. nf90_noerr) THEN
386!              PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
387              PRINT*, "dynetat0: Le champ <"//trim(tnom(iq))//
388     &                "> est absent"
389              PRINT*, "          Il est donc initialise a zero"
390              CALL initial0(ijp1llm,q(1,1,1,iq))
391              nqold=min(iq-1,nqold)
392           ELSE
393           !ierr=nf90_get_var(nid,nvarid,q(1,1,1,iq))
394           ierr=nf90_get_var(nid,nvarid,traceur,corner,edges)
395             IF (ierr .NE. nf90_noerr) THEN
396!                 PRINT*, "dynetat0: Lecture echouee pour "//str3
397               PRINT*, "dynetat0: Lecture echouee pour "//trim(tnom(iq))
398               CALL abort
399             ENDIF
400             q(:,:,:,iq)=traceur(:,:,:)
401           ENDIF
402         ENDDO
403         if ((nqold.lt.nq).and.(nqold.ge.1)) then   
404c        case when new tracer are added in addition to old ones
405             write(*,*)'tracers 1 to ', nqold,'were already present'
406             write(*,*)'tracers ', nqold+1,' to ', nqmx,'are new'
407!             yes=' '
408!            do while ((yes.ne.'y').and.(yes.ne.'n'))
409!             write(*,*) 'Would you like to reindex tracer # 1 ->',nqold
410!             write(*,*) 'to #',nqmx-nqold+1,'->', nqmx,'   (y or n) ?'
411!             read(*,fmt='(a)') yes
412!            end do
413!            if (yes.eq.'y') then
414!              write(*,*) 'OK, let s reindex the tracers'
415!              do l=1,llm
416!                do j=1,jjp1
417!                  do i=1,iip1
418!                    do iq=nqmx,nqmx-nqold+1,-1
419!                       q(i,j,l,iq)=q(i,j,l,iq-nqmx+nqold)   
420!                    end do
421!                    do iq=nqmx-nqold,1,-1
422!                       q(i,j,l,iq)= 0.
423!                    end do
424!                  end do
425!                end do
426!              end do
427!            end if
428         end if
429      ENDIF
430
431! read masse
432      ierr=nf90_inq_varid(nid,"masse",nvarid)
433      IF (ierr .NE. nf90_noerr) THEN
434         PRINT*, "dynetat0: Le champ <masse> est absent"
435         write(*,*)trim(nf90_strerror(ierr))
436         CALL abort
437      ENDIF
438      ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
439      IF (ierr .NE. nf90_noerr) THEN
440         PRINT*, "dynetat0: Lecture echouee pour <masse>"
441         write(*,*)trim(nf90_strerror(ierr))
442         CALL abort
443      ENDIF
444
445! read ps
446      corner(1)=1
447      corner(2)=1
448      corner(3)=indextime
449      edges(1)=iip1
450      edges(2)=jjp1
451      edges(3)=1
452      ierr=nf90_inq_varid(nid,"ps",nvarid)
453      IF (ierr .NE. nf90_noerr) THEN
454         PRINT*, "dynetat0: Le champ <ps> est absent"
455         write(*,*)trim(nf90_strerror(ierr))
456         CALL abort
457      ENDIF
458      ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
459      IF (ierr .NE. nf90_noerr) THEN
460         PRINT*, "dynetat0: Lecture echouee pour <ps>"
461         write(*,*)trim(nf90_strerror(ierr))
462         CALL abort
463      ENDIF
464
465      ierr=nf90_close(nid)
466
467      ! day_ini=day_ini+INT(time) ! obsolete stuff ; 0<time<1 anyways
468      ! time=time-INT(time)
469
470  1   FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem
471     *arrage est differente de la valeur parametree iim =',i4//)
472   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem
473     *arrage est differente de la valeur parametree jjm =',i4//)
474   3  FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema
475     *rrage est differente de la valeur parametree llm =',i4//)
476   4  FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema
477     *rrage est differente de la valeur  dtinteg =',i4//)
478
479      RETURN
480      END
Note: See TracBrowser for help on using the repository browser.