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

Last change on this file since 2306 was 2265, checked in by jyg, 10 years ago

Bug fixes concerning the number "nqo" of water
phases in "q" array: part II (sequel of fixes of
revision 2262).

Modified files:

carbon_cycle_mod.F90
ini_histrac.h
phyredem.F90
phyetat0.F90
traclmdz_mod.F90
write_histrac.h
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.0 KB
Line 
1! $Id: phyetat0.F90 2265 2015-04-18 15:06:36Z jyg $
2
3SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
4
5  USE dimphy, only: klon, zmasq, klev, nslay
6  USE iophy, ONLY : init_iophy_new
7  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
8  USE fonte_neige_mod,  ONLY : fonte_neige_init
9  USE pbl_surface_mod,  ONLY : pbl_surface_init
10  USE surface_data,     ONLY : type_ocean, version_ocean
11  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
12       qsol, fevap, z0m, z0h, agesno, &
13       du_gwd_rando, dv_gwd_rando, entr_therm, f0, fm_therm, &
14       falb_dir, falb_dif, &
15       ftsol, pbl_tke, pctsrf, q_ancien, radpas, radsol, rain_fall, ratqs, &
16       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
17       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
18       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
19       wake_s, zgam, &
20       zmax0, zmea, zpic, zsig, &
21       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
22  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
23  USE infotrac, only: nbtr, nqo, type_trac, tname, niadv
24  USE traclmdz_mod,    ONLY : traclmdz_from_restart
25  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
26  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
27  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
28
29  IMPLICIT none
30  !======================================================================
31  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
32  ! Objet: Lecture de l'etat initial pour la physique
33  !======================================================================
34  include "dimensions.h"
35  include "netcdf.inc"
36  include "dimsoil.h"
37  include "clesphys.h"
38  include "temps.h"
39  include "thermcell.h"
40  include "compbl.h"
41  include "YOMCST.h"
42  !======================================================================
43  CHARACTER*(*) fichnom
44
45  ! les variables globales lues dans le fichier restart
46
47  REAL tsoil(klon, nsoilmx, nbsrf)
48  REAL qsurf(klon, nbsrf)
49  REAL snow(klon, nbsrf)
50  real fder(klon)
51  REAL run_off_lic_0(klon)
52  REAL fractint(klon)
53  REAL trs(klon, nbtr)
54  REAL zts(klon)
55
56  CHARACTER*6 ocean_in
57  LOGICAL ok_veget_in
58
59  INTEGER        longcles
60  PARAMETER    ( longcles = 20 )
61  REAL clesphy0( longcles )
62
63  REAL xmin, xmax
64
65  INTEGER nid, nvarid
66  INTEGER ierr, i, nsrf, isoil , k
67  INTEGER length
68  PARAMETER (length=100)
69  INTEGER it, iiq, isw
70  REAL tab_cntrl(length), tabcntr0(length)
71  CHARACTER*7 str7
72  CHARACTER*2 str2
73  LOGICAL :: found,phyetat0_get,phyetat0_srf
74
75  ! FH1D
76  !     real iolat(jjm+1)
77  real iolat(jjm+1-1/(iim*jjm))
78
79  ! Ouvrir le fichier contenant l'etat initial:
80
81  CALL open_startphy(fichnom)
82
83  ! Lecture des parametres de controle:
84
85  CALL get_var("controle", tab_cntrl)
86
87!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
89  ! Les constantes de la physiques sont lues dans la physique seulement.
90  ! Les egalites du type
91  !             tab_cntrl( 5 )=clesphy0(1)
92  ! sont remplacees par
93  !             clesphy0(1)=tab_cntrl( 5 )
94  ! On inverse aussi la logique.
95  ! On remplit les tab_cntrl avec les parametres lus dans les .def
96!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97
98  DO i = 1, length
99     tabcntr0( i ) = tab_cntrl( i )
100  ENDDO
101
102  tab_cntrl(1)=dtime
103  tab_cntrl(2)=radpas
104
105  ! co2_ppm : value from the previous time step
106  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
107     co2_ppm = tab_cntrl(3)
108     RCO2    = co2_ppm * 1.0e-06  * 44.011/28.97
109     ! ELSE : keep value from .def
110  END IF
111
112  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
113  co2_ppm0   = tab_cntrl(16)
114
115  solaire_etat0      = tab_cntrl(4)
116  tab_cntrl(5)=iflag_con
117  tab_cntrl(6)=nbapp_rad
118
119  if (cycle_diurne) tab_cntrl( 7) =1.
120  if (soil_model) tab_cntrl( 8) =1.
121  if (new_oliq) tab_cntrl( 9) =1.
122  if (ok_orodr) tab_cntrl(10) =1.
123  if (ok_orolf) tab_cntrl(11) =1.
124  if (ok_limitvrai) tab_cntrl(12) =1.
125
126  itau_phy = tab_cntrl(15)
127
128  clesphy0(1)=tab_cntrl( 5 )
129  clesphy0(2)=tab_cntrl( 6 )
130  clesphy0(3)=tab_cntrl( 7 )
131  clesphy0(4)=tab_cntrl( 8 )
132  clesphy0(5)=tab_cntrl( 9 )
133  clesphy0(6)=tab_cntrl( 10 )
134  clesphy0(7)=tab_cntrl( 11 )
135  clesphy0(8)=tab_cntrl( 12 )
136
137  ! Lecture des latitudes (coordonnees):
138
139  CALL get_field("latitude", rlat)
140
141  ! Lecture des longitudes (coordonnees):
142
143  CALL get_field("longitude", rlon)
144
145  ! Lecture du masque terre mer
146
147  CALL get_field("masque", zmasq, found)
148  IF (.NOT. found) THEN
149     PRINT*, 'phyetat0: Le champ <masque> est absent'
150     PRINT *, 'fichier startphy non compatible avec phyetat0'
151  ENDIF
152
153  ! Lecture des fractions pour chaque sous-surface
154
155  ! initialisation des sous-surfaces
156
157  pctsrf = 0.
158
159  ! fraction de terre
160
161  CALL get_field("FTER", pctsrf(:, is_ter), found)
162  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
163
164  ! fraction de glace de terre
165
166  CALL get_field("FLIC", pctsrf(:, is_lic), found)
167  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
168
169  ! fraction d'ocean
170
171  CALL get_field("FOCE", pctsrf(:, is_oce), found)
172  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
173
174  ! fraction glace de mer
175
176  CALL get_field("FSIC", pctsrf(:, is_sic), found)
177  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
178
179  !  Verification de l'adequation entre le masque et les sous-surfaces
180
181  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
182       + pctsrf(1 : klon, is_lic)
183  DO i = 1 , klon
184     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
185        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
186             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
187             , pctsrf(i, is_lic)
188        WRITE(*, *) 'Je force la coherence zmasq=fractint'
189        zmasq(i) = fractint(i)
190     ENDIF
191  END DO
192  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
193       + pctsrf(1 : klon, is_sic)
194  DO i = 1 , klon
195     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
196        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
197             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
198             , pctsrf(i, is_sic)
199        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
200        zmasq(i) = 1. - fractint(i)
201     ENDIF
202  END DO
203
204!===================================================================
205! Lecture des temperatures du sol:
206!===================================================================
207
208  found=phyetat0_get(1,ftsol(:,1),"TS","Surface temperature",283.)
209  IF (found) THEN
210     DO nsrf=2,nbsrf
211        ftsol(:,nsrf)=ftsol(:,1)
212     ENDDO
213  ELSE
214     found=phyetat0_srf(1,ftsol,"TS","Surface temperature",283.)
215  ENDIF
216
217!===================================================================
218  ! Lecture des albedo difus et direct
219!===================================================================
220
221  DO nsrf = 1, nbsrf
222     DO isw=1, nsw
223        IF (isw.GT.99) THEN
224           PRINT*, "Trop de bandes SW"
225           call abort_gcm("phyetat0", "", 1)
226        ENDIF
227        WRITE(str2, '(i2.2)') isw
228        found=phyetat0_srf(1,falb_dir(:, isw,:),"A_dir_SW"//str2//"srf","Direct Albedo",0.2)
229        found=phyetat0_srf(1,falb_dif(:, isw,:),"A_dif_SW"//str2//"srf","Direct Albedo",0.2)
230     ENDDO
231  ENDDO
232
233!===================================================================
234  ! Lecture des temperatures du sol profond:
235!===================================================================
236
237   DO isoil=1, nsoilmx
238        IF (isoil.GT.99) THEN
239           PRINT*, "Trop de couches "
240           call abort_gcm("phyetat0", "", 1)
241        ENDIF
242        WRITE(str2,'(i2.2)') isoil
243        found=phyetat0_srf(1,tsoil(:, isoil,:),"Tsoil"//str2//"srf","Temp soil",0.)
244        IF (.NOT. found) THEN
245           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
246           PRINT*, "          Il prend donc la valeur de surface"
247           tsoil(:, isoil, :)=ftsol(:, :)
248        ENDIF
249   ENDDO
250
251!=======================================================================
252! Lecture precipitation/evaporation
253!=======================================================================
254
255  found=phyetat0_srf(1,qsurf,"QS","Near surface hmidity",0.)
256  found=phyetat0_get(1,qsol,"QSOL","Surface hmidity / bucket",0.)
257  found=phyetat0_srf(1,snow,"SNOW","Surface snow",0.)
258  found=phyetat0_srf(1,fevap,"EVAP","evaporation",0.)
259  found=phyetat0_get(1,snow_fall,"snow_f","snow fall",0.)
260  found=phyetat0_get(1,rain_fall,"rain_f","rain fall",0.)
261
262!=======================================================================
263! Radiation
264!=======================================================================
265
266  found=phyetat0_get(1,solsw,"solsw","net SW radiation surf",0.)
267  found=phyetat0_get(1,sollw,"sollw","net LW radiation surf",0.)
268  found=phyetat0_get(1,sollwdown,"sollwdown","down LW radiation surf",0.)
269  IF (.NOT. found) THEN
270     sollwdown = 0. ;  zts=0.
271     do nsrf=1,nbsrf
272        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
273     enddo
274     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
275  ENDIF
276
277  found=phyetat0_get(1,radsol,"RADS","Solar radiation",0.)
278  found=phyetat0_get(1,fder,"fder","Flux derivative",0.)
279
280
281  ! Lecture de la longueur de rugosite
282  found=phyetat0_srf(1,z0m,"RUG","Z0m ancien",0.001)
283  IF (found) THEN
284     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
285  ELSE
286     found=phyetat0_srf(1,z0m,"Z0m","Roughness length, momentum ",0.001)
287     found=phyetat0_srf(1,z0h,"Z0h","Roughness length, enthalpy ",0.001)
288  ENDIF
289
290  ! Lecture de l'age de la neige:
291  found=phyetat0_srf(1,agesno,"AGESNO","SNOW AGE",0.001)
292
293  ancien_ok=.true.
294  ancien_ok=ancien_ok.AND.phyetat0_get(klev,t_ancien,"TANCIEN","TANCIEN",0.)
295  ancien_ok=ancien_ok.AND.phyetat0_get(klev,q_ancien,"QANCIEN","QANCIEN",0.)
296  ancien_ok=ancien_ok.AND.phyetat0_get(klev,u_ancien,"UANCIEN","UANCIEN",0.)
297  ancien_ok=ancien_ok.AND.phyetat0_get(klev,v_ancien,"VANCIEN","VANCIEN",0.)
298
299  found=phyetat0_get(klev,clwcon,"CLWCON","CLWCON",0.)
300  found=phyetat0_get(klev,rnebcon,"RNEBCON","RNEBCON",0.)
301  found=phyetat0_get(klev,ratqs,"RATQS","RATQS",0.)
302
303  found=phyetat0_get(1,run_off_lic_0,"RUNOFFLIC0","RUNOFFLIC0",0.)
304
305!==================================
306!  TKE
307!==================================
308!
309  IF (iflag_pbl>1) then
310     found=phyetat0_srf(klev+1,pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)
311  ENDIF
312
313  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1 ) then
314    found=phyetat0_srf(klev+1,wake_delta_pbl_tke,"DELTATKE","Del TKE wk/env",0.)
315    found=phyetat0_srf(1,delta_tsurf,"DELTA_TSURF","Delta Ts wk/env ",0.)
316  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )
317
318!==================================
319!  thermiques, poches, convection
320!==================================
321
322! Emanuel
323  found=phyetat0_get(klev,sig1,"sig1","sig1",0.)
324  found=phyetat0_get(klev,w01,"w01","w01",0.)
325
326! Wake
327  found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
328  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
329  found=phyetat0_get(1,wake_s,"WAKE_S","WAKE_S",0.)
330  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
331  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
332  found=phyetat0_get(1,wake_fip,"WAKE_FIP","WAKE_FIP",0.)
333
334! Thermiques
335  found=phyetat0_get(1,zmax0,"ZMAX0","ZMAX0",40.)
336  found=phyetat0_get(1,f0,"F0","F0",1.e-5)
337  found=phyetat0_get(klev+1,fm_therm,"FM_THERM","Thermals mass flux",0.)
338  found=phyetat0_get(klev,entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
339  found=phyetat0_get(klev,detr_therm,"DETR_THERM","Thermals Detrain.",0.)
340
341! ALE/ALP
342  found=phyetat0_get(1,ale_bl,"ALE_BL","ALE BL",0.)
343  found=phyetat0_get(1,ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
344  found=phyetat0_get(1,alp_bl,"ALP_BL","ALP BL",0.)
345
346!===========================================
347  ! Read and send field trs to traclmdz
348!===========================================
349
350  IF (type_trac == 'lmdz') THEN
351     DO it=1, nbtr                                                                 
352!!        iiq=niadv(it+2)                                                           ! jyg
353        iiq=niadv(it+nqo)                                                           ! jyg
354        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
355              "Surf trac"//tname(iiq),0.)
356     END DO
357     CALL traclmdz_from_restart(trs)
358
359     IF (carbon_cycle_cpl) THEN
360        ALLOCATE(co2_send(klon), stat=ierr)
361        IF (ierr /= 0) CALL abort_gcm('phyetat0', 'pb allocation co2_send', 1)
362        found=phyetat0_get(1,co2_send,"co2_send","co2 send",0.)
363     END IF
364  END IF
365
366!===========================================
367!  ondes de gravite / relief
368!===========================================
369
370!  ondes de gravite non orographiques
371  if (ok_gwd_rando) then
372     found=phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
373     found=phyetat0_get(klev,dv_gwd_rando,"dv_gwd_rando","dv_gwd_rando",0.)
374  end if
375
376!  prise en compte du relief sous-maille
377  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
378  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
379  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
380  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
381  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
382  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
383  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
384  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
385  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
386
387!===========================================
388! Initialize ocean
389!===========================================
390
391  IF ( type_ocean == 'slab' ) THEN
392      CALL ocean_slab_init(dtime, pctsrf)
393      found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
394      IF (.NOT. found) THEN
395          PRINT*, "phyetat0: Le champ <tslab> est absent"
396          PRINT*, "Initialisation a tsol_oce"
397          DO i=1,nslay
398              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
399          END DO
400      END IF
401
402      ! Sea ice variables
403      found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
404      IF (version_ocean == 'sicINT') THEN
405          IF (.NOT. found) THEN
406              PRINT*, "phyetat0: Le champ <tice> est absent"
407              PRINT*, "Initialisation a tsol_sic"
408                  tice(:)=ftsol(:,is_sic)
409          END IF
410          IF (.NOT. found) THEN
411              PRINT*, "phyetat0: Le champ <seaice> est absent"
412              PRINT*, "Initialisation a 0/1m suivant fraction glace"
413              seaice(:)=0.
414              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
415                  seaice=917.
416              END WHERE
417          END IF
418      END IF !sea ice INT
419  END IF ! Slab       
420
421  ! on ferme le fichier
422  CALL close_startphy
423
424  ! Initialize module pbl_surface_mod
425
426  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
427
428  ! Initialize module ocean_cpl_mod for the case of coupled ocean
429  IF ( type_ocean == 'couple' ) THEN
430     CALL ocean_cpl_init(dtime, rlon, rlat)
431  ENDIF
432
433  CALL init_iophy_new(rlat, rlon)
434
435  ! Initilialize module fonte_neige_mod     
436  CALL fonte_neige_init(run_off_lic_0)
437
438END SUBROUTINE phyetat0
439
440!===================================================================
441FUNCTION phyetat0_get(nlev,field,name,descr,default)
442!===================================================================
443! Lecture d'un champ avec contrôle
444! Function logique dont le resultat indique si la lecture
445! s'est bien passée
446! On donne une valeur par defaut dans le cas contraire
447!===================================================================
448
449USE iostart, ONLY : get_field
450USE dimphy, only: klon
451
452IMPLICIT NONE
453INCLUDE "iniprint.h"
454
455LOGICAL phyetat0_get
456
457! arguments
458INTEGER,INTENT(IN) :: nlev
459CHARACTER*(*),INTENT(IN) :: name,descr
460REAL,INTENT(IN) :: default
461REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
462
463! Local variables
464LOGICAL found
465
466   CALL get_field(name, field, found)
467   IF (.NOT. found) THEN
468     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
469     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
470     field(:,:)=default
471   ENDIF
472   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
473   phyetat0_get=found
474
475RETURN
476END FUNCTION phyetat0_get
477
478!================================================================
479FUNCTION phyetat0_srf(nlev,field,name,descr,default)
480!===================================================================
481! Lecture d'un champ par sous-surface avec contrôle
482! Function logique dont le resultat indique si la lecture
483! s'est bien passée
484! On donne une valeur par defaut dans le cas contraire
485!===================================================================
486
487USE iostart, ONLY : get_field
488USE dimphy, only: klon
489USE indice_sol_mod, only: nbsrf
490
491IMPLICIT NONE
492INCLUDE "iniprint.h"
493
494LOGICAL phyetat0_srf
495! arguments
496INTEGER,INTENT(IN) :: nlev
497CHARACTER*(*),INTENT(IN) :: name,descr
498REAL,INTENT(IN) :: default
499REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
500
501! Local variables
502LOGICAL found,phyetat0_get
503INTEGER nsrf
504CHARACTER*2 str2
505 
506     IF (nbsrf.GT.99) THEN
507        WRITE(lunout,*) "Trop de sous-mailles"
508        call abort_gcm("phyetat0", "", 1)
509     ENDIF
510
511     DO nsrf = 1, nbsrf
512        WRITE(str2, '(i2.2)') nsrf
513        found= phyetat0_get(nlev,field(:,:, nsrf), &
514        name//str2,descr//" srf:"//str2,default)
515     ENDDO
516
517     phyetat0_srf=found
518
519RETURN
520END FUNCTION phyetat0_srf
521
Note: See TracBrowser for help on using the repository browser.