source: LMDZ5/trunk/libf/phylmd/phyetat0.F90 @ 2052

Last change on this file since 2052 was 2052, checked in by musat, 10 years ago

Correction d'un bogue initialisation sous-surfaces
FC/IM

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