source: LMDZ5/branches/testing/libf/phylmd/phyetat0.F90 @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago

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