source: trunk/LMDZ.COMMON/libf/dyn3dpar/dynetat0.F @ 1243

Last change on this file since 1243 was 1107, checked in by emillour, 11 years ago

Common dynamics: Updates and modifications to enable running Mars physics with

LMDZ.COMMON dynamics:

  • For compilation: adapted makelmdz, create_make_gcm and makelmdz_fcm, bld.cfg to compile aeronomy routines in "aerono$physique" if it exists, and added "-P -traditional" preprocessing flags in "arch-linux-ifort*"
  • Added function "cbrt.F" (cubic root) in 'bibio'
  • Adapted the reading/writing of dynamics (re)start.nc files for Mars. The main issue is that different information (on time, reference and current) is stored and used differently, hence a few if (planet_type =="mars") here and there. Moreover in the martian case there is the possibility to store fields over multiple times. Some Mars-specific variables (ecritphy,ecritstart,timestart) added in control_mod.F and (hour_ini) in temps.h

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