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

Last change on this file since 2252 was 2252, checked in by fhourdin, 9 years ago

Poursuite du nettoyage de phyetat0, et retour à 1+1=2.
Clean and clean, there will always remain soemthing out of it.
  • 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: 16.8 KB
Line 
1! $Id: phyetat0.F90 2252 2015-03-27 16:35:52Z 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 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, 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)
353        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
354              "Surf trac"//tname(iiq),0.)
355     END DO
356     CALL traclmdz_from_restart(trs)
357
358     IF (carbon_cycle_cpl) THEN
359        ALLOCATE(co2_send(klon), stat=ierr)
360        IF (ierr /= 0) CALL abort_gcm('phyetat0', 'pb allocation co2_send', 1)
361        found=phyetat0_get(1,co2_send,"co2_send","co2 send",0.)
362     END IF
363  END IF
364
365!===========================================
366!  ondes de gravite / relief
367!===========================================
368
369!  ondes de gravite non orographiques
370  if (ok_gwd_rando) then
371     found=phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
372     found=phyetat0_get(klev,dv_gwd_rando,"dv_gwd_rando","dv_gwd_rando",0.)
373  end if
374
375!  prise en compte du relief sous-maille
376  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
377  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
378  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
379  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
380  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
381  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
382  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
383  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
384  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
385
386!===========================================
387! Initialize ocean
388!===========================================
389
390  IF ( type_ocean == 'slab' ) THEN
391      CALL ocean_slab_init(dtime, pctsrf)
392      found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
393      IF (.NOT. found) THEN
394          PRINT*, "phyetat0: Le champ <tslab> est absent"
395          PRINT*, "Initialisation a tsol_oce"
396          DO i=1,nslay
397              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
398          END DO
399      END IF
400
401      ! Sea ice variables
402      found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
403      IF (version_ocean == 'sicINT') THEN
404          IF (.NOT. found) THEN
405              PRINT*, "phyetat0: Le champ <tice> est absent"
406              PRINT*, "Initialisation a tsol_sic"
407                  tice(:)=ftsol(:,is_sic)
408          END IF
409          IF (.NOT. found) THEN
410              PRINT*, "phyetat0: Le champ <seaice> est absent"
411              PRINT*, "Initialisation a 0/1m suivant fraction glace"
412              seaice(:)=0.
413              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
414                  seaice=917.
415              END WHERE
416          END IF
417      END IF !sea ice INT
418  END IF ! Slab       
419
420  ! on ferme le fichier
421  CALL close_startphy
422
423  ! Initialize module pbl_surface_mod
424
425  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
426
427  ! Initialize module ocean_cpl_mod for the case of coupled ocean
428  IF ( type_ocean == 'couple' ) THEN
429     CALL ocean_cpl_init(dtime, rlon, rlat)
430  ENDIF
431
432  CALL init_iophy_new(rlat, rlon)
433
434  ! Initilialize module fonte_neige_mod     
435  CALL fonte_neige_init(run_off_lic_0)
436
437END SUBROUTINE phyetat0
438
439!===================================================================
440FUNCTION phyetat0_get(nlev,field,name,descr,default)
441!===================================================================
442! Lecture d'un champ avec contrôle
443! Function logique dont le resultat indique si la lecture
444! s'est bien passée
445! On donne une valeur par defaut dans le cas contraire
446!===================================================================
447
448USE iostart, ONLY : get_field
449USE dimphy, only: klon
450
451IMPLICIT NONE
452INCLUDE "iniprint.h"
453
454LOGICAL phyetat0_get
455
456! arguments
457INTEGER,INTENT(IN) :: nlev
458CHARACTER*(*),INTENT(IN) :: name,descr
459REAL,INTENT(IN) :: default
460REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
461
462! Local variables
463LOGICAL found
464
465   CALL get_field(name, field, found)
466   IF (.NOT. found) THEN
467     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
468     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
469     field(:,:)=default
470   ENDIF
471   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
472   phyetat0_get=found
473
474RETURN
475END FUNCTION phyetat0_get
476
477!================================================================
478FUNCTION phyetat0_srf(nlev,field,name,descr,default)
479!===================================================================
480! Lecture d'un champ par sous-surface avec contrôle
481! Function logique dont le resultat indique si la lecture
482! s'est bien passée
483! On donne une valeur par defaut dans le cas contraire
484!===================================================================
485
486USE iostart, ONLY : get_field
487USE dimphy, only: klon
488USE indice_sol_mod, only: nbsrf
489
490IMPLICIT NONE
491INCLUDE "iniprint.h"
492
493LOGICAL phyetat0_srf
494! arguments
495INTEGER,INTENT(IN) :: nlev
496CHARACTER*(*),INTENT(IN) :: name,descr
497REAL,INTENT(IN) :: default
498REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
499
500! Local variables
501LOGICAL found,phyetat0_get
502INTEGER nsrf
503CHARACTER*2 str2
504 
505     IF (nbsrf.GT.99) THEN
506        WRITE(lunout,*) "Trop de sous-mailles"
507        call abort_gcm("phyetat0", "", 1)
508     ENDIF
509
510     DO nsrf = 1, nbsrf
511        WRITE(str2, '(i2.2)') nsrf
512        found= phyetat0_get(nlev,field(:,:, nsrf), &
513        name//str2,descr//" srf:"//str2,default)
514     ENDDO
515
516     phyetat0_srf=found
517
518RETURN
519END FUNCTION phyetat0_srf
520
Note: See TracBrowser for help on using the repository browser.