source: trunk/LMDZ.MARS/libf/phymars/phyetat0.F @ 190

Last change on this file since 190 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

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