source: LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/phyetat0.F @ 4934

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