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

Last change on this file since 2251 was 2251, checked in by fhourdin, 10 years ago

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