source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/phyetat0.F90 @ 3825

Last change on this file since 3825 was 3825, checked in by ymipsl, 10 years ago

Reorganize geometry and grid modules. Prepare physics for unstructutured grid support. Simplify initialization of physics from dynamic.
Compiled only with dynd3dmem, but not tested for moment.

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