source: trunk/LMDZ.PLUTO.old/libf/phypluto/phyetat0.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

File size: 16.7 KB
Line 
1      SUBROUTINE phyetat0 (fichnom,tab0,Lmodif,nsoil,nq,
2     .           day_ini,time,
3     .           tsurf,tsoil,emis,q2,qsurf)
4      use planet_h         
5      use comgeomfi_h
6      implicit none
7c======================================================================
8c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
9c  Heritage Adaptation to Mars : Yann Wanherdrick
10c Objet: Lecture de l etat initial pour la physique
11c======================================================================
12#include "netcdf.inc"
13#include "dimensions.h"
14#include "dimphys.h"
15#include "surfdat.h"
16#include "comcstfi.h"
17#include"advtrac.h"
18
19c======================================================================
20      INTEGER nbsrf !
21      PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
22!======================================================================
23!  Arguments:
24!  ---------
25!  inputs:
26      character*(*) fichnom ! "startfi.nc" file
27      integer tab0
28      integer Lmodif
29      integer nsoil ! # of soil layers
30      integer ngrid
31      integer nq
32      integer day_ini
33      real time
34
35!  outputs:
36      real tsurf(ngridmx,nbsrf) ! surface temperature
37      real tsoil(ngridmx,nsoil,nbsrf) ! soil temperature
38      real emis(ngridmx) ! surface emissivity
39      real q2(ngridmx, llm+1) !
40      real qsurf(ngridmx,nq) ! tracers on surface
41
42!======================================================================
43!  Local variables:
44
45!      INTEGER radpas
46!      REAL solaire
47
48      real xmin,xmax ! to display min and max of a field
49c
50      INTEGER ig,iq,lmax
51      INTEGER nid, nvarid
52      INTEGER ierr, i, nsrf
53!      integer isoil
54!      INTEGER length
55!      PARAMETER (length=100)
56      CHARACTER*7 str7
57      CHARACTER*2 str2
58      CHARACTER*1 yes
59c
60      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
61      INTEGER nqold
62
63! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
64      logical :: oldtracernames=.false.
65      integer :: count
66      character(len=30) :: txt ! to store some text
67
68 
69c
70c Ouvrir le fichier contenant l etat initial:
71c
72
73      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
74      IF (ierr.NE.NF_NOERR) THEN
75        write(6,*)' Pb d''ouverture du fichier '//fichnom
76        CALL ABORT
77      ENDIF
78
79! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
80!                    qsurf02, ...)
81      count=0
82      do iq=1,nqmx
83        txt= " "
84        write(txt,'(a5,i2.2)')'qsurf',iq
85        ierr=NF_INQ_VARID(nid,txt,nvarid)
86        if (ierr.ne.NF_NOERR) then
87          ! did not find old tracer name
88          exit ! might as well stop here
89        else
90          ! found old tracer name
91          count=count+1
92        endif
93      enddo
94      if (count.eq.nqmx) then
95        write(*,*) "phyetat0:tracers seem to follow old naming ",
96     &             "convention (qsurf01,qsurf02,...)"
97        write(*,*) "   => will work for now ... "
98        write(*,*) "      but you should run newstart to rename them"
99        oldtracernames=.true.
100      endif
101
102c modifications possibles des variables de tab_cntrl
103     
104      PRINT*
105      write(*,*) 'TABFI de phyeta0',Lmodif,tab0
106      call tabfi (nid,Lmodif,tab0,day_ini,lmax,p_rad,
107     .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
108c
109c Lecture des latitudes (coordonnees):
110c
111      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
112      IF (ierr.NE.NF_NOERR) THEN
113         PRINT*, 'phyetat0: Le champ <latitude> est absent'
114         CALL abort
115      ENDIF
116#ifdef NC_DOUBLE
117      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lati)
118#else
119      ierr = NF_GET_VAR_REAL(nid, nvarid, lati)
120#endif
121      IF (ierr.NE.NF_NOERR) THEN
122         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
123         CALL abort
124      ENDIF
125c
126c Lecture des longitudes (coordonnees):
127c
128      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
129      IF (ierr.NE.NF_NOERR) THEN
130         PRINT*, 'phyetat0: Le champ <longitude> est absent'
131         CALL abort
132      ENDIF
133#ifdef NC_DOUBLE
134      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, long)
135#else
136      ierr = NF_GET_VAR_REAL(nid, nvarid, long)
137#endif
138      IF (ierr.NE.NF_NOERR) THEN
139         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
140         CALL abort
141      ENDIF
142c
143c Lecture des aires des mailles:
144c
145      ierr = NF_INQ_VARID (nid, "area", nvarid)
146      IF (ierr.NE.NF_NOERR) THEN
147         PRINT*, 'phyetat0: Le champ <area> est absent'
148         CALL abort
149      ENDIF
150#ifdef NC_DOUBLE
151      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, area)
152#else
153      ierr = NF_GET_VAR_REAL(nid, nvarid, area)
154#endif
155      IF (ierr.NE.NF_NOERR) THEN
156         PRINT*, 'phyetat0: Lecture echouee pour <area>'
157         CALL abort
158      ENDIF
159      xmin = 1.0E+20
160      xmax = -1.0E+20
161      xmin = MINVAL(area)
162      xmax = MAXVAL(area)
163      PRINT*,'Aires des mailles <area>:', xmin, xmax
164c
165c Lecture du geopotentiel au sol:
166c
167      ierr = NF_INQ_VARID (nid, "phisfi", nvarid)
168      IF (ierr.NE.NF_NOERR) THEN
169         PRINT*, 'phyetat0: Le champ <phisfi> est absent'
170         CALL abort
171      ENDIF
172#ifdef NC_DOUBLE
173      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisfi)
174#else
175      ierr = NF_GET_VAR_REAL(nid, nvarid, phisfi)
176#endif
177      IF (ierr.NE.NF_NOERR) THEN
178         PRINT*, 'phyetat0: Lecture echouee pour <phisfi>'
179         CALL abort
180      ENDIF
181      xmin = 1.0E+20
182      xmax = -1.0E+20
183      xmin = MINVAL(phisfi)
184      xmax = MAXVAL(phisfi)
185      PRINT*,'Geopotentiel au sol <phisfi>:', xmin, xmax
186c
187c Lecture de l''albedo du sol nu:
188c
189      ierr = NF_INQ_VARID (nid, "albedodat", nvarid)
190      IF (ierr.NE.NF_NOERR) THEN
191         PRINT*, 'phyetat0: Le champ <albedodat> est absent'
192         CALL abort
193      ENDIF
194#ifdef NC_DOUBLE
195      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albedodat)
196#else
197      ierr = NF_GET_VAR_REAL(nid, nvarid, albedodat)
198#endif
199      IF (ierr.NE.NF_NOERR) THEN
200         PRINT*, 'phyetat0: Lecture echouee pour <albedodat>'
201         CALL abort
202      ENDIF
203      xmin = 1.0E+20
204      xmax = -1.0E+20
205      xmin = MINVAL(albedodat)
206      xmax = MAXVAL(albedodat)
207      PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
208c
209c Lecture de l''inertie thermique du sol:
210c
211!      ierr = NF_INQ_VARID (nid, "inertiedat", nvarid)
212!      IF (ierr.NE.NF_NOERR) THEN
213!         PRINT*, 'phyetat0: Le champ <inertiedat> est absent'
214!         CALL abort
215!      ENDIF
216!#ifdef NC_DOUBLE
217!      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, inertiedat)
218!#else
219!      ierr = NF_GET_VAR_REAL(nid, nvarid, inertiedat)
220!#endif
221!      IF (ierr.NE.NF_NOERR) THEN
222!         PRINT*, 'phyetat0: Lecture echouee pour <inertiedat>'
223!         CALL abort
224!      ENDIF
225!      xmin = 1.0E+20
226!      xmax = -1.0E+20
227!      xmin = MINVAL(inertiedat)
228!      xmax = MAXVAL(inertiedat)
229!      PRINT*,'Inertie thermique du sol <inertiedat>:', xmin, xmax
230c
231c ZMEA
232c
233      ierr = NF_INQ_VARID (nid, "ZMEA", nvarid)
234      IF (ierr.NE.NF_NOERR) THEN
235         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
236         CALL abort
237      ENDIF
238#ifdef NC_DOUBLE
239      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmea)
240#else
241      ierr = NF_GET_VAR_REAL(nid, nvarid, zmea)
242#endif
243      IF (ierr.NE.NF_NOERR) THEN
244         PRINT*, 'phyetat0: Lecture echouee pour <ZMEA>'
245         CALL abort
246      ENDIF
247      xmin = 1.0E+20
248      xmax = -1.0E+20
249      DO i = 1, ngridmx
250         xmin = MIN(zmea(i),xmin)
251         xmax = MAX(zmea(i),xmax)
252      ENDDO
253      PRINT*,'<zmea>:', xmin, xmax
254c
255c ZSTD
256c
257      ierr = NF_INQ_VARID (nid, "ZSTD", nvarid)
258      IF (ierr.NE.NF_NOERR) THEN
259         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
260         CALL abort
261      ENDIF
262#ifdef NC_DOUBLE
263      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zstd)
264#else
265      ierr = NF_GET_VAR_REAL(nid, nvarid, zstd)
266#endif
267      IF (ierr.NE.NF_NOERR) THEN
268         PRINT*, 'phyetat0: Lecture echouee pour <ZSTD>'
269         CALL abort
270      ENDIF
271      xmin = 1.0E+20
272      xmax = -1.0E+20
273      DO i = 1, ngridmx
274         xmin = MIN(zstd(i),xmin)
275         xmax = MAX(zstd(i),xmax)
276      ENDDO
277      PRINT*,'<zstd>:', xmin, xmax
278c
279c ZSIG
280c
281      ierr = NF_INQ_VARID (nid, "ZSIG", nvarid)
282      IF (ierr.NE.NF_NOERR) THEN
283         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
284         CALL abort
285      ENDIF
286#ifdef NC_DOUBLE
287      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zsig)
288#else
289      ierr = NF_GET_VAR_REAL(nid, nvarid, zsig)
290#endif
291      IF (ierr.NE.NF_NOERR) THEN
292         PRINT*, 'phyetat0: Lecture echouee pour <ZSIG>'
293         CALL abort
294      ENDIF
295      xmin = 1.0E+20
296      xmax = -1.0E+20
297      DO i = 1, ngridmx
298         xmin = MIN(zsig(i),xmin)
299         xmax = MAX(zsig(i),xmax)
300      ENDDO
301      PRINT*,'<zsig>:', xmin, xmax
302c
303c ZGAM
304c
305      ierr = NF_INQ_VARID (nid, "ZGAM", nvarid)
306      IF (ierr.NE.NF_NOERR) THEN
307         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
308         CALL abort
309      ENDIF
310#ifdef NC_DOUBLE
311      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zgam)
312#else
313      ierr = NF_GET_VAR_REAL(nid, nvarid, zgam)
314#endif
315      IF (ierr.NE.NF_NOERR) THEN
316         PRINT*, 'phyetat0: Lecture echouee pour <ZGAM>'
317         CALL abort
318      ENDIF
319      xmin = 1.0E+20
320      xmax = -1.0E+20
321      DO i = 1, ngridmx
322         xmin = MIN(zgam(i),xmin)
323         xmax = MAX(zgam(i),xmax)
324      ENDDO
325      PRINT*,'<zgam>:', xmin, xmax
326c
327c ZTHE
328c
329      ierr = NF_INQ_VARID (nid, "ZTHE", nvarid)
330      IF (ierr.NE.NF_NOERR) THEN
331         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
332         CALL abort
333      ENDIF
334#ifdef NC_DOUBLE
335      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zthe)
336#else
337      ierr = NF_GET_VAR_REAL(nid, nvarid, zthe)
338#endif
339      IF (ierr.NE.NF_NOERR) THEN
340         PRINT*, 'phyetat0: Lecture echouee pour <ZTHE>'
341         CALL abort
342      ENDIF
343      xmin = 1.0E+20
344      xmax = -1.0E+20
345      DO i = 1, ngridmx
346         xmin = MIN(zthe(i),xmin)
347         xmax = MAX(zthe(i),xmax)
348      ENDDO
349      PRINT*,'<zthe>:', xmin, xmax
350
351c
352c Lecture des temperatures du sol:
353c
354      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
355      IF (ierr.NE.NF_NOERR) THEN
356         PRINT*, 'phyetat0: Le champ <tsurf> est absent'
357         PRINT*, '          Mais je vais essayer de lire TS**'
358         IF (nbsrf.GT.99) THEN
359            PRINT*, "Trop de sous-mailles"
360            CALL abort
361         ENDIF
362         DO nsrf = 1, nbsrf
363           WRITE(str2,'(i2.2)') nsrf
364           ierr = NF_INQ_VARID (nid, "TS"//str2, nvarid)
365           IF (ierr.NE.NF_NOERR) THEN
366           PRINT*, "phyetat0: Le champ <TS"//str2//"> est absent"
367              CALL abort
368           ENDIF
369#ifdef NC_DOUBLE
370           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsurf(1,nsrf))
371#else
372           ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(1,nsrf))
373#endif
374           IF (ierr.NE.NF_NOERR) THEN
375             PRINT*, "phyetat0: Lecture echouee pour <TS"//str2//">"
376             CALL abort
377           ENDIF
378           xmin = 1.0E+20
379           xmax = -1.0E+20
380           xmin = MINVAL(tsurf)
381           xmax = MAXVAL(tsurf)
382           PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
383         ENDDO
384      ELSE
385         PRINT*, 'phyetat0: Le champ <tsurf> est present'
386         PRINT*, '          J ignore donc les autres temperatures TS**'
387#ifdef NC_DOUBLE
388         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsurf(1,1))
389#else
390         ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(1,1))
391#endif
392         IF (ierr.NE.NF_NOERR) THEN
393            PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
394            CALL abort
395         ENDIF
396         xmin = 1.0E+20
397         xmax = -1.0E+20
398         xmin = MINVAL(tsurf)
399         xmax = MAXVAL(tsurf)
400         PRINT*,'Temperature du sol <tsurf>', xmin, xmax
401         IF (nbsrf >= 2) THEN
402            DO nsrf = 2, nbsrf
403               DO i = 1, ngridmx
404                  tsurf(i,nsrf) = tsurf(i,1)
405               ENDDO
406            ENDDO
407         ENDIF
408      ENDIF
409c
410c Lecture des temperatures du sol profond:
411c
412!      IF (nsoil.GT.99 .OR. nbsrf.GT.99) THEN
413!         PRINT*, "Trop de couches ou sous-mailles"
414!         CALL abort
415!      ENDIF
416!      DO nsrf = 1, nbsrf
417!         DO isoil=1, nsoil
418!            WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
419!            ierr = NF_INQ_VARID (nid, 'tsoil', nvarid)
420!            IF (ierr.NE.NF_NOERR) THEN
421!               PRINT*, "phyetat0: Le champ <tsoil> est absent"
422!               PRINT*, "          Il prend donc la valeur de surface"
423!               DO i=1, ngridmx
424!                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
425!               ENDDO
426!            ELSE
427!#ifdef NC_DOUBLE
428!              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,1,nsrf))
429!#else
430!              ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,1,nsrf))
431!#endif
432!               IF (ierr.NE.NF_NOERR) THEN
433!                  PRINT*, "Lecture echouee pour <tsoil>"
434!                  CALL abort
435!               ENDIF
436!            ENDIF
437!         ENDDO
438!      ENDDO
439!      xmin = 1.0E+20
440!      xmax = -1.0E+20
441!      xmin = MINVAL(tsoil)
442!      xmax = MAXVAL(tsoil)
443!      PRINT*,'Temperatures du sol profond <tsoil>', xmin, xmax
444c
445c Surface emissivity
446c
447      ierr = NF_INQ_VARID (nid, "emis", nvarid)
448      IF (ierr.NE.NF_NOERR) THEN
449         PRINT*, 'phyetat0: Le champ <emis> est absent'
450         CALL abort
451      ENDIF
452#ifdef NC_DOUBLE
453      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emis)
454#else
455      ierr = NF_GET_VAR_REAL(nid, nvarid, emis)
456#endif
457      IF (ierr.NE.NF_NOERR) THEN
458         PRINT*, 'phyetat0: Lecture echouee pour <emis>'
459         CALL abort
460      ENDIF
461      xmin = 1.0E+20
462      xmax = -1.0E+20
463      xmin = MINVAL(emis)
464      xmax = MAXVAL(emis)
465      PRINT*,'Min/Max Surface emissivity <emis>:', xmin, xmax
466c
467c pbl wind variance
468c
469      ierr = NF_INQ_VARID (nid, "q2", nvarid)
470      IF (ierr.NE.NF_NOERR) THEN
471         PRINT*, 'phyetat0: Le champ <q2> est absent'
472         CALL abort
473      ENDIF
474#ifdef NC_DOUBLE
475      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q2)
476#else
477      ierr = NF_GET_VAR_REAL(nid, nvarid, q2)
478#endif
479      IF (ierr.NE.NF_NOERR) THEN
480         PRINT*, 'phyetat0: Lecture echouee pour <q2>'
481         CALL abort
482      ENDIF
483      xmin = 1.0E+20
484      xmax = -1.0E+20
485      xmin = MINVAL(q2)
486      xmax = MAXVAL(q2)
487      PRINT*,'pbl wind variance <q2>:', xmin, xmax
488c
489c tracer on surface
490c
491
492      IF(nq.GE.1) THEN
493         nqold=nq
494         DO iq=1,nq
495!            str7(1:5)='qsurf'
496!            WRITE(str7(6:7),'(i2.2)') iq
497!            ierr = NF_INQ_VARID (nid,str7,nvarid)
498           IF (oldtracernames) THEN
499             txt=" "
500             write(txt,'(a5,i2.2)')'qsurf',iq
501           ELSE
502             txt=tnom(iq)
503             if (txt.eq."h2o_vap") then
504               ! There is no surface tracer for h2o_vap;
505               ! "h2o_ice" should be loaded instead
506               txt="h2o_ice"
507               write(*,*) 'phyetat0: loading surface tracer',
508     &                     ' h2o_ice instead of h2o_vap'
509             endif
510           ENDIF ! of IF (oldtracernames) THEN
511           ierr=NF_INQ_VARID(nid,txt,nvarid)
512           IF (ierr.NE.NF_NOERR) THEN
513             write(*,*) 'PHYETAT0: WARNING : surface tracer',trim(txt),
514     &                  ' not found in file'
515             write(*,*) trim(txt), ' set to 0'
516             do ig=1,ngridmx
517               qsurf(ig,iq)=0.
518             end do
519             nqold=min(iq-1,nqold)
520           ELSE
521#ifdef NC_DOUBLE
522             ierr = NF_GET_VAR_DOUBLE(nid, nvarid,qsurf(1,iq))
523#else
524             ierr = NF_GET_VAR_REAL(nid, nvarid,qsurf(1,iq))
525#endif
526             IF (ierr.NE.NF_NOERR) THEN
527               PRINT*, 'phyetat0: Lecture echouee pour <',trim(txt),'>'
528               CALL abort
529             ENDIF
530           ENDIF
531           xmin = 1.0E+20
532           xmax = -1.0E+20
533           xmin = MINVAL(qsurf(1:ngridmx,iq))
534           xmax = MAXVAL(qsurf(1:ngridmx,iq))
535           PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
536         ENDDO
537         if ((nqold.lt.nq).and.(nqold.ge.1)) then
538c        case when new tracer are added in addition to old ones
539             write(*,*)'qsurf 1 to ', nqold,'were already present'
540             write(*,*)'qsurf ', nqold+1,' to ', nqmx,'are new'
541             write(*,*)' and initialized to zero'
542             qsurf(:,nqold+1:nqmx)=0.0
543!            yes=' '
544!            do while ((yes.ne.'y').and.(yes.ne.'n'))
545!             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
546!             write(*,*) 'to #',nqmx-nqold+1,'->', nqmx,'   (y or n) ?'
547!             read(*,fmt='(a)') yes
548!            end do
549!            if (yes.eq.'y') then
550!              write(*,*) 'OK, let s reindex qsurf'
551!                 do ig=1,ngridmx
552!                    do iq=nqmx,nqmx-nqold+1,-1
553!                       qsurf(ig,iq)=qsurf(ig,iq-nqmx+nqold)
554!                    end do
555!                    do iq=nqmx-nqold,1,-1
556!                       qsurf(ig,iq)= 0.
557!                    end do
558!                 end do
559!            end if
560         end if ! of if ((nqold.lt.nq).and.(nqold.ge.1))
561      ENDIF ! of IF(nq.GE.1)
562
563! Call to soil_settings, in order to read soil temperatures,
564! as well as thermal inertia and volumetric heat capacity
565
566      call soil_settings(nid,ngridmx,nsoil,tsurf,tsoil)
567c
568c Fermer le fichier:
569c
570      ierr = NF_CLOSE(nid)
571c
572      RETURN
573      END
Note: See TracBrowser for help on using the repository browser.