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

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

Temporary disabled printing min and max value for orography parametrization, causing floating point exception.

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