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

Last change on this file since 3838 was 3838, checked in by millour, 9 years ago

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