source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phyetat0.F90 @ 3547

Last change on this file since 3547 was 3411, checked in by Laurent Fairhead, 6 years ago

Undoing merge with trunk (r3356) to properly register Yann's latest modifications
  • 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 Id
File size: 21.1 KB
Line 
1! $Id: phyetat0.F90 3411 2018-11-05 14:24:59Z jbmadeleine $
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, phys_tstep, &
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  ELSE
338    ! initialize treedrg(), because it will be written in restartphy.nc
339    treedrg(:,:,:) = 0.0
340  ENDIF
341
342  ! Lecture de l'age de la neige:
343  found=phyetat0_srf(1,agesno,"AGESNO","SNOW AGE",0.001)
344
345  ancien_ok=.true.
346  ancien_ok=ancien_ok.AND.phyetat0_get(klev,t_ancien,"TANCIEN","TANCIEN",0.)
347  ancien_ok=ancien_ok.AND.phyetat0_get(klev,q_ancien,"QANCIEN","QANCIEN",0.)
348  ancien_ok=ancien_ok.AND.phyetat0_get(klev,ql_ancien,"QLANCIEN","QLANCIEN",0.)
349  ancien_ok=ancien_ok.AND.phyetat0_get(klev,qs_ancien,"QSANCIEN","QSANCIEN",0.)
350  ancien_ok=ancien_ok.AND.phyetat0_get(klev,u_ancien,"UANCIEN","UANCIEN",0.)
351  ancien_ok=ancien_ok.AND.phyetat0_get(klev,v_ancien,"VANCIEN","VANCIEN",0.)
352  ancien_ok=ancien_ok.AND.phyetat0_get(1,prw_ancien,"PRWANCIEN","PRWANCIEN",0.)
353  ancien_ok=ancien_ok.AND.phyetat0_get(1,prlw_ancien,"PRLWANCIEN","PRLWANCIEN",0.)
354  ancien_ok=ancien_ok.AND.phyetat0_get(1,prsw_ancien,"PRSWANCIEN","PRSWANCIEN",0.)
355
356  ! Ehouarn: addtional tests to check if t_ancien, q_ancien contain
357  !          dummy values (as is the case when generated by ce0l,
358  !          or by iniaqua)
359  if ( (maxval(q_ancien).eq.minval(q_ancien))       .or. &
360       (maxval(ql_ancien).eq.minval(ql_ancien))     .or. &
361       (maxval(qs_ancien).eq.minval(qs_ancien))     .or. &
362       (maxval(prw_ancien).eq.minval(prw_ancien))   .or. &
363       (maxval(prlw_ancien).eq.minval(prlw_ancien)) .or. &
364       (maxval(prsw_ancien).eq.minval(prsw_ancien)) .or. &
365       (maxval(t_ancien).eq.minval(t_ancien)) ) then
366    ancien_ok=.false.
367  endif
368
369  found=phyetat0_get(klev,clwcon,"CLWCON","CLWCON",0.)
370  found=phyetat0_get(klev,rnebcon,"RNEBCON","RNEBCON",0.)
371  found=phyetat0_get(klev,ratqs,"RATQS","RATQS",0.)
372
373  found=phyetat0_get(1,run_off_lic_0,"RUNOFFLIC0","RUNOFFLIC0",0.)
374
375!==================================
376!  TKE
377!==================================
378!
379  IF (iflag_pbl>1) then
380     found=phyetat0_srf(klev+1,pbl_tke,"TKE","Turb. Kinetic. Energ. ",1.e-8)
381  ENDIF
382
383  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1 ) then
384    found=phyetat0_srf(klev+1,wake_delta_pbl_tke,"DELTATKE","Del TKE wk/env",0.)
385    found=phyetat0_srf(1,delta_tsurf,"DELTA_TSURF","Delta Ts wk/env ",0.)
386  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )
387
388!==================================
389!  thermiques, poches, convection
390!==================================
391
392! Emanuel
393  found=phyetat0_get(klev,sig1,"sig1","sig1",0.)
394  found=phyetat0_get(klev,w01,"w01","w01",0.)
395
396! Wake
397  found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
398  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
399  found=phyetat0_get(1,wake_s,"WAKE_S","Wake frac. area",0.)
400!jyg<
401!  Set wake_dens to -1000. when there is no restart so that the actual
402!  initialization is made in calwake.
403!!  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
404  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",-1000.)
405!>jyg
406  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
407  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
408  found=phyetat0_get(1,wake_fip,"WAKE_FIP","WAKE_FIP",0.)
409
410! Thermiques
411  found=phyetat0_get(1,zmax0,"ZMAX0","ZMAX0",40.)
412  found=phyetat0_get(1,f0,"F0","F0",1.e-5)
413  found=phyetat0_get(klev+1,fm_therm,"FM_THERM","Thermals mass flux",0.)
414  found=phyetat0_get(klev,entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
415  found=phyetat0_get(klev,detr_therm,"DETR_THERM","Thermals Detrain.",0.)
416
417! ALE/ALP
418  found=phyetat0_get(1,ale_bl,"ALE_BL","ALE BL",0.)
419  found=phyetat0_get(1,ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
420  found=phyetat0_get(1,alp_bl,"ALP_BL","ALP BL",0.)
421
422!===========================================
423  ! Read and send field trs to traclmdz
424!===========================================
425
426  IF (type_trac == 'lmdz') THEN
427     DO it=1, nbtr                                                                 
428!!        iiq=niadv(it+2)                                                           ! jyg
429        iiq=niadv(it+nqo)                                                           ! jyg
430        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
431              "Surf trac"//tname(iiq),0.)
432     END DO
433     CALL traclmdz_from_restart(trs)
434
435     IF (carbon_cycle_cpl) THEN
436        ALLOCATE(co2_send(klon), stat=ierr)
437        IF (ierr /= 0) CALL abort_physic('phyetat0', 'pb allocation co2_send', 1)
438        found=phyetat0_get(1,co2_send,"co2_send","co2 send",0.)
439     END IF
440  END IF
441
442!===========================================
443!  ondes de gravite / relief
444!===========================================
445
446!  ondes de gravite non orographiques
447  if (ok_gwd_rando) found = &
448       phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
449  IF (.not. ok_hines .and. ok_gwd_rando) found &
450       = phyetat0_get(klev,du_gwd_front,"du_gwd_front","du_gwd_front",0.)
451
452!  prise en compte du relief sous-maille
453  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
454  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
455  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
456  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
457  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
458  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
459  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
460  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
461  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
462
463!===========================================
464! Initialize ocean
465!===========================================
466
467  IF ( type_ocean == 'slab' ) THEN
468      CALL ocean_slab_init(phys_tstep, pctsrf)
469      IF (nslay.EQ.1) THEN
470        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
471        IF (.NOT. found) THEN
472            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
473        END IF
474      ELSE
475          DO i=1,nslay
476            WRITE(str2,'(i2.2)') i
477            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.) 
478          END DO
479      END IF
480      IF (.NOT. found) THEN
481          PRINT*, "phyetat0: Le champ <tslab> est absent"
482          PRINT*, "Initialisation a tsol_oce"
483          DO i=1,nslay
484              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
485          END DO
486      END IF
487
488      ! Sea ice variables
489      IF (version_ocean == 'sicINT') THEN
490          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
491          IF (.NOT. found) THEN
492              PRINT*, "phyetat0: Le champ <tice> est absent"
493              PRINT*, "Initialisation a tsol_sic"
494                  tice(:)=ftsol(:,is_sic)
495          END IF
496          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
497          IF (.NOT. found) THEN
498              PRINT*, "phyetat0: Le champ <seaice> est absent"
499              PRINT*, "Initialisation a 0/1m suivant fraction glace"
500              seaice(:)=0.
501              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
502                  seaice=917.
503              END WHERE
504          END IF
505      END IF !sea ice INT
506  END IF ! Slab       
507
508  ! on ferme le fichier
509  CALL close_startphy
510
511  ! Initialize module pbl_surface_mod
512
513  CALL pbl_surface_init(fder, snow, qsurf, tsoil)
514
515  ! Initialize module ocean_cpl_mod for the case of coupled ocean
516  IF ( type_ocean == 'couple' ) THEN
517     CALL ocean_cpl_init(phys_tstep, longitude_deg, latitude_deg)
518  ENDIF
519
520!  CALL init_iophy_new(latitude_deg, longitude_deg)
521
522  ! Initilialize module fonte_neige_mod     
523  CALL fonte_neige_init(run_off_lic_0)
524
525END SUBROUTINE phyetat0
526
527!===================================================================
528FUNCTION phyetat0_get(nlev,field,name,descr,default)
529!===================================================================
530! Lecture d'un champ avec contrôle
531! Function logique dont le resultat indique si la lecture
532! s'est bien passée
533! On donne une valeur par defaut dans le cas contraire
534!===================================================================
535
536USE iostart, ONLY : get_field
537USE dimphy, only: klon
538USE print_control_mod, ONLY: lunout
539
540IMPLICIT NONE
541
542LOGICAL phyetat0_get
543
544! arguments
545INTEGER,INTENT(IN) :: nlev
546CHARACTER*(*),INTENT(IN) :: name,descr
547REAL,INTENT(IN) :: default
548REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
549
550! Local variables
551LOGICAL found
552
553   CALL get_field(name, field, found)
554   IF (.NOT. found) THEN
555     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
556     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
557     field(:,:)=default
558   ENDIF
559   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
560   phyetat0_get=found
561
562RETURN
563END FUNCTION phyetat0_get
564
565!================================================================
566FUNCTION phyetat0_srf(nlev,field,name,descr,default)
567!===================================================================
568! Lecture d'un champ par sous-surface avec contrôle
569! Function logique dont le resultat indique si la lecture
570! s'est bien passée
571! On donne une valeur par defaut dans le cas contraire
572!===================================================================
573
574USE iostart, ONLY : get_field
575USE dimphy, only: klon
576USE indice_sol_mod, only: nbsrf
577USE print_control_mod, ONLY: lunout
578
579IMPLICIT NONE
580
581LOGICAL phyetat0_srf
582! arguments
583INTEGER,INTENT(IN) :: nlev
584CHARACTER*(*),INTENT(IN) :: name,descr
585REAL,INTENT(IN) :: default
586REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
587
588! Local variables
589LOGICAL found,phyetat0_get
590INTEGER nsrf
591CHARACTER*2 str2
592 
593     IF (nbsrf.GT.99) THEN
594        WRITE(lunout,*) "Trop de sous-mailles"
595        call abort_physic("phyetat0", "", 1)
596     ENDIF
597
598     DO nsrf = 1, nbsrf
599        WRITE(str2, '(i2.2)') nsrf
600        found= phyetat0_get(nlev,field(:,:, nsrf), &
601        name//str2,descr//" srf:"//str2,default)
602     ENDDO
603
604     phyetat0_srf=found
605
606RETURN
607END FUNCTION phyetat0_srf
608
Note: See TracBrowser for help on using the repository browser.