source: LMDZ5/branches/LMDZ_tree_FC/libf/phylmd/phyetat0.F90 @ 5371

Last change on this file since 5371 was 2940, checked in by fcheruy, 7 years ago

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