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

Last change on this file since 3983 was 3960, checked in by ymipsl, 8 years ago

Fix uninitialized variables in LMDZ physics.

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