source: trunk/MESOSCALE/LMDZ.MARS/libf_gcm/phymars/phyetat0.F @ 1242

Last change on this file since 1242 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

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