source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/phyetat0.F @ 3536

Last change on this file since 3536 was 1319, checked in by Laurent Fairhead, 15 years ago
  • Modifications to the start and limit creation routines to account for different

calendars

  • Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

  • Force readaerosol.F90 to read in aerosols file with 12 timesteps

  • Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

  • Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

  • Forcer readaerosol à ne lire que des fichiers à 12 pas de temps
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 28.9 KB
Line 
1!
2! $Header$
3!
4c
5c
6      SUBROUTINE phyetat0 (fichnom,
7     .           clesphy0,
8     .           tabcntr0)
9
10      USE dimphy
11      USE mod_grid_phy_lmdz
12      USE mod_phys_lmdz_para
13      USE iophy
14      USE ocean_cpl_mod,    ONLY : ocean_cpl_init
15      USE fonte_neige_mod,  ONLY : fonte_neige_init
16      USE pbl_surface_mod,  ONLY : pbl_surface_init
17      USE surface_data,     ONLY : type_ocean
18      USE phys_state_var_mod
19      USE iostart
20      USE write_field_phy
21      USE infotrac
22      USE traclmdz_mod,    ONLY : traclmdz_from_restart
23      USE carbon_cycle_mod,ONLY : carbon_cycle_tr, carbon_cycle_cpl
24
25      IMPLICIT none
26c======================================================================
27c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
28c Objet: Lecture de l'etat initial pour la physique
29c======================================================================
30#include "dimensions.h"
31#include "netcdf.inc"
32#include "indicesol.h"
33#include "dimsoil.h"
34#include "clesphys.h"
35#include "temps.h"
36#include "thermcell.h"
37#include "compbl.h"
38c======================================================================
39      CHARACTER*(*) fichnom
40
41c les variables globales lues dans le fichier restart
42
43      REAL tsoil(klon,nsoilmx,nbsrf)
44      REAL tslab(klon), seaice(klon)
45      REAL qsurf(klon,nbsrf)
46      REAL qsol(klon)
47      REAL snow(klon,nbsrf)
48      REAL evap(klon,nbsrf)
49      real fder(klon)
50      REAL frugs(klon,nbsrf)
51      REAL agesno(klon,nbsrf)
52      REAL run_off_lic_0(klon)
53      REAL fractint(klon)
54      REAL trs(klon,nbtr)
55
56      CHARACTER*6 ocean_in
57      LOGICAL ok_veget_in
58
59      INTEGER        longcles
60      PARAMETER    ( longcles = 20 )
61      REAL clesphy0( longcles )
62c
63      REAL xmin, xmax
64c
65      INTEGER nid, nvarid
66      INTEGER ierr, i, nsrf, isoil ,k
67      INTEGER length
68      PARAMETER (length=100)
69      INTEGER it, iiq
70      REAL tab_cntrl(length), tabcntr0(length)
71      CHARACTER*7 str7
72      CHARACTER*2 str2
73      LOGICAL :: found
74
75c FH1D
76c     real iolat(jjm+1)
77      real iolat(jjm+1-1/iim)
78c
79c Ouvrir le fichier contenant l'etat initial:
80c
81
82     
83      CALL open_startphy(fichnom)
84     
85
86c
87c Lecture des parametres de controle:
88c
89      CALL get_var("controle",tab_cntrl)
90       
91c
92!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
93! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
94! Les constantes de la physiques sont lues dans la physique seulement.
95! Les egalites du type
96!             tab_cntrl( 5 )=clesphy0(1)
97! sont remplacees par
98!             clesphy0(1)=tab_cntrl( 5 )
99! On inverse aussi la logique.
100! On remplit les tab_cntrl avec les parametres lus dans les .def
101!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
102
103         DO i = 1, length
104           tabcntr0( i ) = tab_cntrl( i )
105         ENDDO
106c
107         tab_cntrl(1)=dtime
108         tab_cntrl(2)=radpas
109
110c co2_ppm : value from the previous time step
111         IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
112            co2_ppm = tab_cntrl(3)
113            RCO2    = co2_ppm * 1.0e-06  * 44.011/28.97
114c ELSE : keep value from .def
115         END IF
116
117c co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
118         co2_ppm0   = tab_cntrl(16)
119
120         solaire_etat0      = tab_cntrl(4)
121         tab_cntrl(5)=iflag_con
122         tab_cntrl(6)=nbapp_rad
123
124         if (cycle_diurne) tab_cntrl( 7) =1.
125         if (soil_model) tab_cntrl( 8) =1.
126         if (new_oliq) tab_cntrl( 9) =1.
127         if (ok_orodr) tab_cntrl(10) =1.
128         if (ok_orolf) tab_cntrl(11) =1.
129         if (ok_limitvrai) tab_cntrl(12) =1.
130
131
132      itau_phy = tab_cntrl(15)
133
134       
135
136         IF( clesphy0(1).NE.tab_cntrl( 5 ) )  THEN
137             clesphy0(1)=tab_cntrl( 5 )
138         ENDIF
139
140         IF( clesphy0(2).NE.tab_cntrl( 6 ) )  THEN
141             clesphy0(2)=tab_cntrl( 6 )
142         ENDIF
143
144         IF( clesphy0(3).NE.tab_cntrl( 7 ) )  THEN
145             clesphy0(3)=tab_cntrl( 7 )
146         ENDIF
147
148         IF( clesphy0(4).NE.tab_cntrl( 8 ) )  THEN
149             clesphy0(4)=tab_cntrl( 8 )
150         ENDIF
151
152         IF( clesphy0(5).NE.tab_cntrl( 9 ) )  THEN
153             clesphy0(5)=tab_cntrl( 9 )
154         ENDIF
155
156         IF( clesphy0(6).NE.tab_cntrl( 10 ) )  THEN
157             clesphy0(6)=tab_cntrl( 10 )
158         ENDIF
159
160         IF( clesphy0(7).NE.tab_cntrl( 11 ) )  THEN
161             clesphy0(7)=tab_cntrl( 11 )
162         ENDIF
163
164         IF( clesphy0(8).NE.tab_cntrl( 12 ) )  THEN
165             clesphy0(8)=tab_cntrl( 12 )
166         ENDIF
167
168
169c
170c Lecture des latitudes (coordonnees):
171c
172      CALL get_field("latitude",rlat)
173
174c
175c Lecture des longitudes (coordonnees):
176c
177      CALL get_field("longitude",rlon)
178
179C
180C
181C Lecture du masque terre mer
182C
183      CALL get_field("masque",zmasq,found)
184      IF (.NOT. found) THEN
185        PRINT*, 'phyetat0: Le champ <masque> est absent'
186        PRINT *, 'fichier startphy non compatible avec phyetat0'
187      ENDIF
188
189       
190C Lecture des fractions pour chaque sous-surface
191C
192C initialisation des sous-surfaces
193C
194      pctsrf = 0.
195C
196C fraction de terre
197C
198
199      CALL get_field("FTER",pctsrf(:,is_ter),found)
200      IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
201
202C
203C fraction de glace de terre
204C
205      CALL get_field("FLIC",pctsrf(:,is_lic),found)
206      IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
207
208C
209C fraction d'ocean
210C
211      CALL get_field("FOCE",pctsrf(:,is_oce),found)
212      IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
213
214C
215C fraction glace de mer
216C
217      CALL get_field("FSIC",pctsrf(:,is_sic),found)
218      IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
219
220C
221C  Verification de l'adequation entre le masque et les sous-surfaces
222C
223      fractint( 1 : klon) = pctsrf(1 : klon, is_ter)
224     $    + pctsrf(1 : klon, is_lic)
225      DO i = 1 , klon
226        IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
227            WRITE(*,*) 'phyetat0: attention fraction terre pas ',
228     $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
229     $          ,pctsrf(i, is_lic)
230            zmasq(i) = fractint(i)
231        ENDIF
232      END DO
233      fractint (1 : klon) =  pctsrf(1 : klon, is_oce)
234     $    + pctsrf(1 : klon, is_sic)
235      DO i = 1 , klon
236        IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
237            WRITE(*,*) 'phyetat0 attention fraction ocean pas ',
238     $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
239     $          ,pctsrf(i, is_sic)
240            zmasq(i) = fractint(i)
241        ENDIF
242      END DO
243
244C
245c Lecture des temperatures du sol:
246c
247
248       CALL get_field("TS",ftsol(:,1),found)
249       IF (.NOT. found) THEN
250         PRINT*, 'phyetat0: Le champ <TS> est absent'
251         PRINT*, '          Mais je vais essayer de lire TS**'
252         DO nsrf = 1, nbsrf
253           IF (nsrf.GT.99) THEN
254             PRINT*, "Trop de sous-mailles"
255             CALL abort
256           ENDIF
257           WRITE(str2,'(i2.2)') nsrf
258           CALL get_field("TS"//str2,ftsol(:,nsrf))
259
260           xmin = 1.0E+20
261           xmax = -1.0E+20
262           DO i = 1, klon
263              xmin = MIN(ftsol(i,nsrf),xmin)
264              xmax = MAX(ftsol(i,nsrf),xmax)
265           ENDDO
266           PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
267         ENDDO
268      ELSE
269         PRINT*, 'phyetat0: Le champ <TS> est present'
270         PRINT*, '          J ignore donc les autres temperatures TS**'
271         xmin = 1.0E+20
272         xmax = -1.0E+20
273         DO i = 1, klon
274            xmin = MIN(ftsol(i,1),xmin)
275            xmax = MAX(ftsol(i,1),xmax)
276         ENDDO
277         PRINT*,'Temperature du sol <TS>', xmin, xmax
278         DO nsrf = 2, nbsrf
279         DO i = 1, klon
280            ftsol(i,nsrf) = ftsol(i,1)
281         ENDDO
282         ENDDO
283      ENDIF
284
285c
286c Lecture des temperatures du sol profond:
287c
288      DO nsrf = 1, nbsrf
289        DO isoil=1, nsoilmx
290          IF (isoil.GT.99 .AND. nsrf.GT.99) THEN
291            PRINT*, "Trop de couches ou sous-mailles"
292            CALL abort
293          ENDIF
294          WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
295         
296          CALL get_field('Tsoil'//str7,tsoil(:,isoil,nsrf),found)
297          IF (.NOT. found) THEN
298            PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
299            PRINT*, "          Il prend donc la valeur de surface"
300            DO i=1, klon
301               tsoil(i,isoil,nsrf)=ftsol(i,nsrf)
302            ENDDO
303          ENDIF
304        ENDDO
305      ENDDO
306c
307c Lecture de l'humidite de l'air juste au dessus du sol:
308c
309
310      CALL get_field("QS",qsurf(:,1),found)
311      IF (.NOT. found) THEN
312         PRINT*, 'phyetat0: Le champ <QS> est absent'
313         PRINT*, '          Mais je vais essayer de lire QS**'
314         DO nsrf = 1, nbsrf
315           IF (nsrf.GT.99) THEN
316             PRINT*, "Trop de sous-mailles"
317             CALL abort
318           ENDIF
319           WRITE(str2,'(i2.2)') nsrf
320           CALL get_field("QS"//str2,qsurf(:,nsrf))
321           xmin = 1.0E+20
322           xmax = -1.0E+20
323           DO i = 1, klon
324              xmin = MIN(qsurf(i,nsrf),xmin)
325              xmax = MAX(qsurf(i,nsrf),xmax)
326           ENDDO
327           PRINT*,'Humidite pres du sol QS**:', nsrf, xmin, xmax
328         ENDDO
329      ELSE
330         PRINT*, 'phyetat0: Le champ <QS> est present'
331         PRINT*, '          J ignore donc les autres humidites QS**'
332         xmin = 1.0E+20
333         xmax = -1.0E+20
334         DO i = 1, klon
335            xmin = MIN(qsurf(i,1),xmin)
336            xmax = MAX(qsurf(i,1),xmax)
337         ENDDO
338         PRINT*,'Humidite pres du sol <QS>', xmin, xmax
339         DO nsrf = 2, nbsrf
340           DO i = 1, klon
341             qsurf(i,nsrf) = qsurf(i,1)
342           ENDDO
343         ENDDO
344      ENDIF
345
346C
347C Eau dans le sol (pour le modele de sol "bucket")
348C
349      CALL get_field("QSOL",qsol,found)
350      IF (.NOT. found) THEN
351        PRINT*, 'phyetat0: Le champ <QSOL> est absent'
352        PRINT*, '          Valeur par defaut nulle'
353          qsol(:)=0.
354      ENDIF
355
356      xmin = 1.0E+20
357      xmax = -1.0E+20
358      DO i = 1, klon
359        xmin = MIN(qsol(i),xmin)
360        xmax = MAX(qsol(i),xmax)
361      ENDDO
362      PRINT*,'Eau dans le sol (mm) <QSOL>', xmin, xmax
363
364c
365c Lecture de neige au sol:
366c
367
368      CALL get_field("SNOW",snow(:,1),found)
369      IF (.NOT. found) THEN
370        PRINT*, 'phyetat0: Le champ <SNOW> est absent'
371        PRINT*, '          Mais je vais essayer de lire SNOW**'
372        DO nsrf = 1, nbsrf
373          IF (nsrf.GT.99) THEN
374            PRINT*, "Trop de sous-mailles"
375            CALL abort
376          ENDIF
377          WRITE(str2,'(i2.2)') nsrf
378          CALL get_field( "SNOW"//str2,snow(:,nsrf))
379          xmin = 1.0E+20
380          xmax = -1.0E+20
381          DO i = 1, klon
382            xmin = MIN(snow(i,nsrf),xmin)
383            xmax = MAX(snow(i,nsrf),xmax)
384          ENDDO
385          PRINT*,'Neige du sol SNOW**:', nsrf, xmin, xmax
386        ENDDO
387      ELSE
388         PRINT*, 'phyetat0: Le champ <SNOW> est present'
389         PRINT*, '          J ignore donc les autres neiges SNOW**'
390         xmin = 1.0E+20
391         xmax = -1.0E+20
392         DO i = 1, klon
393            xmin = MIN(snow(i,1),xmin)
394            xmax = MAX(snow(i,1),xmax)
395         ENDDO
396         PRINT*,'Neige du sol <SNOW>', xmin, xmax
397         DO nsrf = 2, nbsrf
398         DO i = 1, klon
399            snow(i,nsrf) = snow(i,1)
400         ENDDO
401         ENDDO
402      ENDIF
403c
404c Lecture de albedo de l'interval visible au sol:
405c
406      CALL get_field("ALBE",falb1(:,1),found)
407      IF (.NOT. found) THEN
408         PRINT*, 'phyetat0: Le champ <ALBE> est absent'
409         PRINT*, '          Mais je vais essayer de lire ALBE**'
410         DO nsrf = 1, nbsrf
411           IF (nsrf.GT.99) THEN
412             PRINT*, "Trop de sous-mailles"
413             CALL abort
414           ENDIF
415           WRITE(str2,'(i2.2)') nsrf
416           CALL get_field("ALBE"//str2,falb1(:,nsrf))
417           xmin = 1.0E+20
418           xmax = -1.0E+20
419           DO i = 1, klon
420              xmin = MIN(falb1(i,nsrf),xmin)
421              xmax = MAX(falb1(i,nsrf),xmax)
422           ENDDO
423           PRINT*,'Albedo du sol ALBE**:', nsrf, xmin, xmax
424         ENDDO
425      ELSE
426         PRINT*, 'phyetat0: Le champ <ALBE> est present'
427         PRINT*, '          J ignore donc les autres ALBE**'
428         xmin = 1.0E+20
429         xmax = -1.0E+20
430         DO i = 1, klon
431            xmin = MIN(falb1(i,1),xmin)
432            xmax = MAX(falb1(i,1),xmax)
433         ENDDO
434         PRINT*,'Neige du sol <ALBE>', xmin, xmax
435         DO nsrf = 2, nbsrf
436           DO i = 1, klon
437            falb1(i,nsrf) = falb1(i,1)
438           ENDDO
439         ENDDO
440      ENDIF
441
442c
443c Lecture de albedo au sol dans l'interval proche infra-rouge:
444c
445      CALL get_field("ALBLW",falb2(:,1),found)
446      IF (.NOT. found) THEN
447         PRINT*, 'phyetat0: Le champ <ALBLW> est absent'
448         PRINT*, '          Mais je vais prendre ALBE**'
449         DO nsrf = 1, nbsrf
450           DO i = 1, klon
451             falb2(i,nsrf) = falb1(i,nsrf)
452           ENDDO
453         ENDDO
454      ELSE
455         PRINT*, 'phyetat0: Le champ <ALBLW> est present'
456         PRINT*, '          J ignore donc les autres ALBLW**'
457         xmin = 1.0E+20
458         xmax = -1.0E+20
459         DO i = 1, klon
460            xmin = MIN(falb2(i,1),xmin)
461            xmax = MAX(falb2(i,1),xmax)
462         ENDDO
463         PRINT*,'Neige du sol <ALBLW>', xmin, xmax
464         DO nsrf = 2, nbsrf
465           DO i = 1, klon
466             falb2(i,nsrf) = falb2(i,1)
467           ENDDO
468         ENDDO
469      ENDIF
470c
471c Lecture de evaporation: 
472c
473      CALL get_field("EVAP",evap(:,1),found)
474      IF (.NOT. found) THEN
475         PRINT*, 'phyetat0: Le champ <EVAP> est absent'
476         PRINT*, '          Mais je vais essayer de lire EVAP**'
477         DO nsrf = 1, nbsrf
478           IF (nsrf.GT.99) THEN
479             PRINT*, "Trop de sous-mailles"
480             CALL abort
481           ENDIF
482           WRITE(str2,'(i2.2)') nsrf
483           CALL get_field("EVAP"//str2, evap(:,nsrf))
484           xmin = 1.0E+20
485           xmax = -1.0E+20
486           DO i = 1, klon
487              xmin = MIN(evap(i,nsrf),xmin)
488              xmax = MAX(evap(i,nsrf),xmax)
489           ENDDO
490           PRINT*,'evap du sol EVAP**:', nsrf, xmin, xmax
491         ENDDO
492      ELSE
493         PRINT*, 'phyetat0: Le champ <EVAP> est present'
494         PRINT*, '          J ignore donc les autres EVAP**'
495         xmin = 1.0E+20
496         xmax = -1.0E+20
497         DO i = 1, klon
498            xmin = MIN(evap(i,1),xmin)
499            xmax = MAX(evap(i,1),xmax)
500         ENDDO
501         PRINT*,'Evap du sol <EVAP>', xmin, xmax
502         DO nsrf = 2, nbsrf
503         DO i = 1, klon
504            evap(i,nsrf) = evap(i,1)
505         ENDDO
506         ENDDO
507      ENDIF
508c
509c Lecture precipitation liquide:
510c
511      CALL get_field("rain_f",rain_fall)
512      xmin = 1.0E+20
513      xmax = -1.0E+20
514      DO i = 1, klon
515         xmin = MIN(rain_fall(i),xmin)
516         xmax = MAX(rain_fall(i),xmax)
517      ENDDO
518      PRINT*,'Precipitation liquide rain_f:', xmin, xmax
519c
520c Lecture precipitation solide:
521c
522      CALL get_field("snow_f",snow_fall)
523      xmin = 1.0E+20
524      xmax = -1.0E+20
525      DO i = 1, klon
526         xmin = MIN(snow_fall(i),xmin)
527         xmax = MAX(snow_fall(i),xmax)
528      ENDDO
529      PRINT*,'Precipitation solide snow_f:', xmin, xmax
530c
531c Lecture rayonnement solaire au sol:
532c
533      CALL get_field("solsw",solsw,found)
534      IF (.NOT. found) THEN
535         PRINT*, 'phyetat0: Le champ <solsw> est absent'
536         PRINT*, 'mis a zero'
537         solsw(:) = 0.
538      ENDIF
539      xmin = 1.0E+20
540      xmax = -1.0E+20
541      DO i = 1, klon
542         xmin = MIN(solsw(i),xmin)
543         xmax = MAX(solsw(i),xmax)
544      ENDDO
545      PRINT*,'Rayonnement solaire au sol solsw:', xmin, xmax
546c
547c Lecture rayonnement IF au sol:
548c
549      CALL get_field("sollw",sollw,found)
550      IF (.NOT. found) THEN
551         PRINT*, 'phyetat0: Le champ <sollw> est absent'
552         PRINT*, 'mis a zero'
553         sollw = 0.
554      ENDIF
555      xmin = 1.0E+20
556      xmax = -1.0E+20
557      DO i = 1, klon
558         xmin = MIN(sollw(i),xmin)
559         xmax = MAX(sollw(i),xmax)
560      ENDDO
561      PRINT*,'Rayonnement IF au sol sollw:', xmin, xmax
562     
563c
564c Lecture derive des flux:
565c
566      CALL get_field("fder",fder,found)
567      IF (.NOT. found) THEN
568         PRINT*, 'phyetat0: Le champ <fder> est absent'
569         PRINT*, 'mis a zero'
570         fder = 0.
571      ENDIF
572      xmin = 1.0E+20
573      xmax = -1.0E+20
574      DO i = 1, klon
575         xmin = MIN(fder(i),xmin)
576         xmax = MAX(fder(i),xmax)
577      ENDDO
578      PRINT*,'Derive des flux fder:', xmin, xmax
579
580c
581c Lecture du rayonnement net au sol:
582c
583      CALL get_field("RADS",radsol)
584      xmin = 1.0E+20
585      xmax = -1.0E+20
586      DO i = 1, klon
587         xmin = MIN(radsol(i),xmin)
588         xmax = MAX(radsol(i),xmax)
589      ENDDO
590      PRINT*,'Rayonnement net au sol radsol:', xmin, xmax
591c
592c Lecture de la longueur de rugosite
593c
594c
595      CALL get_field("RUG",frugs(:,1),found)
596      IF (.NOT. found) THEN
597         PRINT*, 'phyetat0: Le champ <RUG> est absent'
598         PRINT*, '          Mais je vais essayer de lire RUG**'
599         DO nsrf = 1, nbsrf
600           IF (nsrf.GT.99) THEN
601             PRINT*, "Trop de sous-mailles"
602             CALL abort
603           ENDIF
604           WRITE(str2,'(i2.2)') nsrf
605           CALL get_field("RUG"//str2,frugs(:,nsrf))
606           xmin = 1.0E+20
607           xmax = -1.0E+20
608           DO i = 1, klon
609              xmin = MIN(frugs(i,nsrf),xmin)
610              xmax = MAX(frugs(i,nsrf),xmax)
611           ENDDO
612           PRINT*,'rugosite du sol RUG**:', nsrf, xmin, xmax
613         ENDDO
614      ELSE
615         PRINT*, 'phyetat0: Le champ <RUG> est present'
616         PRINT*, '          J ignore donc les autres RUG**'
617         xmin = 1.0E+20
618         xmax = -1.0E+20
619         DO i = 1, klon
620            xmin = MIN(frugs(i,1),xmin)
621            xmax = MAX(frugs(i,1),xmax)
622         ENDDO
623         PRINT*,'rugosite <RUG>', xmin, xmax
624         DO nsrf = 2, nbsrf
625         DO i = 1, klon
626            frugs(i,nsrf) = frugs(i,1)
627         ENDDO
628         ENDDO
629      ENDIF
630
631c
632c Lecture de l'age de la neige:
633c
634      CALL get_field("AGESNO",agesno(:,1),found)
635      IF (.NOT. found) THEN
636         PRINT*, 'phyetat0: Le champ <AGESNO> est absent'
637         PRINT*, '          Mais je vais essayer de lire AGESNO**'
638         DO nsrf = 1, nbsrf
639           IF (nsrf.GT.99) THEN
640             PRINT*, "Trop de sous-mailles"
641             CALL abort
642           ENDIF
643           WRITE(str2,'(i2.2)') nsrf
644           CALL get_field("AGESNO"//str2,agesno(:,nsrf),found)
645           IF (.NOT. found) THEN
646              PRINT*, "phyetat0: Le champ <AGESNO"//str2//"> est absent"
647              agesno = 50.0
648           ENDIF
649           xmin = 1.0E+20
650           xmax = -1.0E+20
651           DO i = 1, klon
652              xmin = MIN(agesno(i,nsrf),xmin)
653              xmax = MAX(agesno(i,nsrf),xmax)
654           ENDDO
655           PRINT*,'Age de la neige AGESNO**:', nsrf, xmin, xmax
656         ENDDO
657      ELSE
658         PRINT*, 'phyetat0: Le champ <AGESNO> est present'
659         PRINT*, '          J ignore donc les autres AGESNO**'
660         xmin = 1.0E+20
661         xmax = -1.0E+20
662         DO i = 1, klon
663            xmin = MIN(agesno(i,1),xmin)
664            xmax = MAX(agesno(i,1),xmax)
665         ENDDO
666         PRINT*,'Age de la neige <AGESNO>', xmin, xmax
667         DO nsrf = 2, nbsrf
668         DO i = 1, klon
669            agesno(i,nsrf) = agesno(i,1)
670         ENDDO
671         ENDDO
672      ENDIF
673
674c
675      CALL get_field("ZMEA", zmea)
676      xmin = 1.0E+20
677      xmax = -1.0E+20
678      DO i = 1, klon
679         xmin = MIN(zmea(i),xmin)
680         xmax = MAX(zmea(i),xmax)
681      ENDDO
682      PRINT*,'OROGRAPHIE SOUS-MAILLE zmea:', xmin, xmax
683c
684c
685      CALL get_field("ZSTD",zstd)
686      xmin = 1.0E+20
687      xmax = -1.0E+20
688      DO i = 1, klon
689         xmin = MIN(zstd(i),xmin)
690         xmax = MAX(zstd(i),xmax)
691      ENDDO
692      PRINT*,'OROGRAPHIE SOUS-MAILLE zstd:', xmin, xmax
693c
694c
695      CALL get_field("ZSIG",zsig)
696      xmin = 1.0E+20
697      xmax = -1.0E+20
698      DO i = 1, klon
699         xmin = MIN(zsig(i),xmin)
700         xmax = MAX(zsig(i),xmax)
701      ENDDO
702      PRINT*,'OROGRAPHIE SOUS-MAILLE zsig:', xmin, xmax
703c
704c
705      CALL get_field("ZGAM",zgam)
706      xmin = 1.0E+20
707      xmax = -1.0E+20
708      DO i = 1, klon
709         xmin = MIN(zgam(i),xmin)
710         xmax = MAX(zgam(i),xmax)
711      ENDDO
712      PRINT*,'OROGRAPHIE SOUS-MAILLE zgam:', xmin, xmax
713c
714c
715      CALL get_field("ZTHE",zthe)
716      xmin = 1.0E+20
717      xmax = -1.0E+20
718      DO i = 1, klon
719         xmin = MIN(zthe(i),xmin)
720         xmax = MAX(zthe(i),xmax)
721      ENDDO
722      PRINT*,'OROGRAPHIE SOUS-MAILLE zthe:', xmin, xmax
723c
724c
725      CALL get_field("ZPIC",zpic)
726      xmin = 1.0E+20
727      xmax = -1.0E+20
728      DO i = 1, klon
729         xmin = MIN(zpic(i),xmin)
730         xmax = MAX(zpic(i),xmax)
731      ENDDO
732      PRINT*,'OROGRAPHIE SOUS-MAILLE zpic:', xmin, xmax
733c
734      CALL get_field("ZVAL",zval)
735      xmin = 1.0E+20
736      xmax = -1.0E+20
737      DO i = 1, klon
738         xmin = MIN(zval(i),xmin)
739         xmax = MAX(zval(i),xmax)
740      ENDDO
741      PRINT*,'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax
742c
743c
744      CALL get_field("RUGSREL",rugoro)
745      xmin = 1.0E+20
746      xmax = -1.0E+20
747      DO i = 1, klon
748         xmin = MIN(rugoro(i),xmin)
749         xmax = MAX(rugoro(i),xmax)
750      ENDDO
751      PRINT*,'Rugosite relief (ecart-type) rugsrel:', xmin, xmax
752c
753c
754     
755c
756      ancien_ok = .TRUE.
757
758      CALL get_field("TANCIEN",t_ancien,found)
759      IF (.NOT. found) THEN
760         PRINT*, "phyetat0: Le champ <TANCIEN> est absent"
761         PRINT*, "Depart legerement fausse. Mais je continue"
762         ancien_ok = .FALSE.
763      ENDIF
764
765
766      CALL get_field("QANCIEN",q_ancien,found)
767      IF (.NOT. found) THEN
768         PRINT*, "phyetat0: Le champ <QANCIEN> est absent"
769         PRINT*, "Depart legerement fausse. Mais je continue"
770         ancien_ok = .FALSE.
771      ENDIF
772
773      u_ancien = 0.0   !AXC: We don't have u_ancien and v_ancien in the start
774      v_ancien = 0.0   !AXC: files, therefore they have to be initialized.
775c
776
777      clwcon=0.
778      CALL get_field("CLWCON",clwcon(:,1),found)
779      IF (.NOT. found) THEN
780         PRINT*, "phyetat0: Le champ CLWCON est absent"
781         PRINT*, "Depart legerement fausse. Mais je continue"
782      ENDIF
783      xmin = 1.0E+20
784      xmax = -1.0E+20
785      xmin = MINval(clwcon)
786      xmax = MAXval(clwcon)
787      PRINT*,'Eau liquide convective (ecart-type) clwcon:', xmin, xmax
788c
789      rnebcon = 0.
790      CALL get_field("RNEBCON",rnebcon(:,1),found)
791      IF (.NOT. found) THEN
792         PRINT*, "phyetat0: Le champ RNEBCON est absent"
793         PRINT*, "Depart legerement fausse. Mais je continue"
794      ENDIF
795      xmin = 1.0E+20
796      xmax = -1.0E+20
797      xmin = MINval(rnebcon)
798      xmax = MAXval(rnebcon)
799      PRINT*,'Nebulosite convective (ecart-type) rnebcon:', xmin, xmax
800
801c
802c Lecture ratqs
803c
804      ratqs=0.
805      CALL get_field("RATQS",ratqs(:,1),found)
806      IF (.NOT. found) THEN
807         PRINT*, "phyetat0: Le champ <RATQS> est absent"
808         PRINT*, "Depart legerement fausse. Mais je continue"
809      ENDIF
810      xmin = 1.0E+20
811      xmax = -1.0E+20
812      xmin = MINval(ratqs)
813      xmax = MAXval(ratqs)
814      PRINT*,'(ecart-type) ratqs:', xmin, xmax
815c
816c Lecture run_off_lic_0
817c
818      CALL get_field("RUNOFFLIC0",run_off_lic_0,found)
819      IF (.NOT. found) THEN
820         PRINT*, "phyetat0: Le champ <RUNOFFLIC0> est absent"
821         PRINT*, "Depart legerement fausse. Mais je continue"
822         run_off_lic_0 = 0.
823      ENDIF
824      xmin = 1.0E+20
825      xmax = -1.0E+20
826      xmin = MINval(run_off_lic_0)
827      xmax = MAXval(run_off_lic_0)
828      PRINT*,'(ecart-type) run_off_lic_0:', xmin, xmax
829
830
831c Lecture de l'energie cinetique turbulente
832c
833
834      IF (iflag_pbl>1) then
835        DO nsrf = 1, nbsrf
836          IF (nsrf.GT.99) THEN
837            PRINT*, "Trop de sous-mailles"
838            CALL abort
839          ENDIF
840          WRITE(str2,'(i2.2)') nsrf
841          CALL get_field("TKE"//str2,pbl_tke(:,1:klev,nsrf),found)
842          IF (.NOT. found) THEN
843            PRINT*, "phyetat0: <TKE"//str2//"> est absent"
844            pbl_tke(:,:,nsrf)=1.e-8
845          ENDIF
846          xmin = 1.0E+20
847          xmax = -1.0E+20
848          DO k = 1, klev
849            DO i = 1, klon
850              xmin = MIN(pbl_tke(i,k,nsrf),xmin)
851              xmax = MAX(pbl_tke(i,k,nsrf),xmax)
852            ENDDO
853          ENDDO
854          PRINT*,'Temperature du sol TKE**:', nsrf, xmin, xmax
855        ENDDO
856      ENDIF
857c
858c zmax0
859      CALL get_field("ZMAX0",zmax0,found)
860      IF (.NOT. found) THEN
861        PRINT*, "phyetat0: Le champ <ZMAX0> est absent"
862        PRINT*, "Depart legerement fausse. Mais je continue"
863        zmax0=40.
864      ENDIF
865      xmin = 1.0E+20
866      xmax = -1.0E+20
867      xmin = MINval(zmax0)
868      xmax = MAXval(zmax0)
869      PRINT*,'(ecart-type) zmax0:', xmin, xmax
870c
871c           f0(ig)=1.e-5
872c f0
873      CALL get_field("F0",f0,found)
874      IF (.NOT. found) THEN
875         PRINT*, "phyetat0: Le champ <f0> est absent"
876         PRINT*, "Depart legerement fausse. Mais je continue"
877         f0=1.e-5
878      ENDIF
879      xmin = 1.0E+20
880      xmax = -1.0E+20
881      xmin = MINval(f0)
882      xmax = MAXval(f0)
883      PRINT*,'(ecart-type) f0:', xmin, xmax
884c
885c ema_work1
886c
887      CALL get_field("EMA_WORK1",ema_work1,found)
888      IF (.NOT. found) THEN
889        PRINT*, "phyetat0: Le champ <EMA_WORK1> est absent"
890        PRINT*, "Depart legerement fausse. Mais je continue"
891        ema_work1=0.
892      ELSE
893        xmin = 1.0E+20
894        xmax = -1.0E+20
895        DO k = 1, klev
896          DO i = 1, klon
897            xmin = MIN(ema_work1(i,k),xmin)
898            xmax = MAX(ema_work1(i,k),xmax)
899          ENDDO
900        ENDDO
901        PRINT*,'ema_work1:', xmin, xmax
902      ENDIF
903c
904c ema_work2
905c
906      CALL get_field("EMA_WORK2",ema_work2,found)
907      IF (.NOT. found) THEN
908        PRINT*, "phyetat0: Le champ <EMA_WORK2> est absent"
909        PRINT*, "Depart legerement fausse. Mais je continue"
910        ema_work2=0.
911      ELSE
912        xmin = 1.0E+20
913        xmax = -1.0E+20
914        DO k = 1, klev
915          DO i = 1, klon
916            xmin = MIN(ema_work2(i,k),xmin)
917            xmax = MAX(ema_work2(i,k),xmax)
918          ENDDO
919        ENDDO
920        PRINT*,'ema_work2:', xmin, xmax
921      ENDIF
922c
923c wake_deltat
924c
925      CALL get_field("WAKE_DELTAT",wake_deltat,found)
926      IF (.NOT. found) THEN
927        PRINT*, "phyetat0: Le champ <WAKE_DELTAT> est absent"
928        PRINT*, "Depart legerement fausse. Mais je continue"
929        wake_deltat=0.
930      ELSE
931        xmin = 1.0E+20
932        xmax = -1.0E+20
933        DO k = 1, klev
934          DO i = 1, klon
935            xmin = MIN(wake_deltat(i,k),xmin)
936            xmax = MAX(wake_deltat(i,k),xmax)
937          ENDDO
938        ENDDO
939        PRINT*,'wake_deltat:', xmin, xmax
940      ENDIF
941c
942c wake_deltaq
943c   
944      CALL get_field("WAKE_DELTAQ",wake_deltaq,found)
945      IF (.NOT. found) THEN
946        PRINT*, "phyetat0: Le champ <WAKE_DELTAQ> est absent"
947        PRINT*, "Depart legerement fausse. Mais je continue"
948        wake_deltaq=0.
949      ELSE
950        xmin = 1.0E+20
951        xmax = -1.0E+20
952        DO k = 1, klev
953          DO i = 1, klon
954            xmin = MIN(wake_deltaq(i,k),xmin)
955            xmax = MAX(wake_deltaq(i,k),xmax)
956          ENDDO
957        ENDDO
958        PRINT*,'wake_deltaq:', xmin, xmax
959      ENDIF
960c
961c wake_s
962c
963      CALL get_field("WAKE_S",wake_s,found)
964      IF (.NOT. found) THEN
965        PRINT*, "phyetat0: Le champ <WAKE_S> est absent"
966        PRINT*, "Depart legerement fausse. Mais je continue"
967        wake_s=0.
968      ENDIF
969      xmin = 1.0E+20
970      xmax = -1.0E+20
971      xmin = MINval(wake_s)
972      xmax = MAXval(wake_s)
973      PRINT*,'(ecart-type) wake_s:', xmin, xmax
974c
975c wake_cstar
976c
977      CALL get_field("WAKE_CSTAR",wake_cstar,found)
978      IF (.NOT. found) THEN
979         PRINT*, "phyetat0: Le champ <WAKE_CSTAR> est absent"
980         PRINT*, "Depart legerement fausse. Mais je continue"
981         wake_cstar=0.
982      ENDIF
983      xmin = 1.0E+20
984      xmax = -1.0E+20
985      xmin = MINval(wake_cstar)
986      xmax = MAXval(wake_cstar)
987      PRINT*,'(ecart-type) wake_cstar:', xmin, xmax
988c
989c wake_fip
990c
991      CALL get_field("WAKE_FIP",wake_fip,found)
992      IF (.NOT. found) THEN
993         PRINT*, "phyetat0: Le champ <WAKE_FIP> est absent"
994         PRINT*, "Depart legerement fausse. Mais je continue"
995         wake_fip=0.
996      ENDIF
997      xmin = 1.0E+20
998      xmax = -1.0E+20
999      xmin = MINval(wake_fip)
1000      xmax = MAXval(wake_fip)
1001      PRINT*,'(ecart-type) wake_fip:', xmin, xmax
1002c
1003c Read and send field trs to traclmdz
1004c
1005      IF (type_trac == 'lmdz') THEN
1006         DO it=1,nbtr
1007            iiq=niadv(it+2)
1008            CALL get_field("trs_"//tname(iiq),trs(:,it),found)
1009            IF (.NOT. found) THEN
1010               PRINT*,
1011     $           "phyetat0: Le champ <trs_"//tname(iiq)//"> est absent"
1012               PRINT*, "Depart legerement fausse. Mais je continue"
1013               trs(:,it) = 0.
1014            ENDIF
1015            xmin = 1.0E+20
1016            xmax = -1.0E+20
1017            xmin = MINval(trs(:,it))
1018            xmax = MAXval(trs(:,it))
1019            PRINT*,"(ecart-type) trs_"//tname(iiq)//" :", xmin, xmax
1020
1021         END DO
1022         
1023         CALL traclmdz_from_restart(trs)
1024      END IF
1025
1026
1027c on ferme le fichier
1028      CALL close_startphy
1029
1030      CALL init_iophy_new(rlat,rlon)
1031       
1032
1033c
1034c Initialize module pbl_surface_mod
1035c
1036      CALL pbl_surface_init(qsol, fder, snow, qsurf,
1037     $     evap, frugs, agesno, tsoil)
1038
1039c Initialize module ocean_cpl_mod for the case of coupled ocean
1040      IF ( type_ocean == 'couple' ) THEN
1041         CALL ocean_cpl_init(dtime, rlon, rlat)
1042      ENDIF
1043c
1044c Initilialize module fonte_neige_mod     
1045c
1046      CALL fonte_neige_init(run_off_lic_0)
1047
1048
1049      RETURN
1050      END
Note: See TracBrowser for help on using the repository browser.