source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/phyetat0.F90 @ 3773

Last change on this file since 3773 was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 21.5 KB
RevLine 
[3331]1! $Id: phyetat0.F90 2569 2016-06-14 07:16:45Z 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, 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, zgam, zmax0, zmea, zpic, zsig, &
20#ifdef ISO
21       fxtevap, xtsol, xt_ancien, xtl_ancien, xts_ancien, wake_deltaxt, &
22       xtrain_fall,xtsnow_fall, &
23#endif
24       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m
25  USE geometry_mod, ONLY : longitude_deg, latitude_deg
26  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
27  USE infotrac_phy, only: nbtr, nqo, type_trac, tname, niadv
28  USE traclmdz_mod,    ONLY : traclmdz_from_restart
29  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
30  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
31  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
32  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
33#ifdef ISO
34  USE infotrac_phy, ONLY: ntraciso,niso,iso_num,itr_indice
35  USE isotopes_routines_mod, ONLY: phyisoetat0
36  USE isotopes_mod, ONLY: iso_eau
37#ifdef ISOVERIF
38  USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_egalite
39#endif
40#endif
41
42  IMPLICIT none
43  !======================================================================
44  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
45  ! Objet: Lecture de l'etat initial pour la physique
46  !======================================================================
47  include "netcdf.inc"
48  include "dimsoil.h"
49  include "clesphys.h"
50  include "thermcell.h"
51  include "compbl.h"
52  include "YOMCST.h"
53  !======================================================================
54  CHARACTER*(*) fichnom
55
56  ! les variables globales lues dans le fichier restart
57
58  REAL tsoil(klon, nsoilmx, nbsrf)
59  REAL qsurf(klon, nbsrf)
60  REAL snow(klon, nbsrf)
61  real fder(klon)
62  REAL run_off_lic_0(klon)
63  REAL fractint(klon)
64  REAL trs(klon, nbtr)
65  REAL zts(klon)
66
67  CHARACTER*6 ocean_in
68  LOGICAL ok_veget_in
69
70  INTEGER        longcles
71  PARAMETER    ( longcles = 20 )
72  REAL clesphy0( longcles )
73
74  REAL xmin, xmax
75
76  INTEGER nid, nvarid
77  INTEGER ierr, i, nsrf, isoil , k
78  INTEGER length
79  PARAMETER (length=100)
80  INTEGER it, iiq, isw
81  REAL tab_cntrl(length), tabcntr0(length)
82  CHARACTER*7 str7
83  CHARACTER*2 str2
84  LOGICAL :: found,phyetat0_get,phyetat0_srf
85  REAL :: lon_startphy(klon), lat_startphy(klon)
86
87#ifdef ISO
88  REAL xtsnow(niso,klon, nbsrf)
89  REAL xtrun_off_lic_0(niso,klon)
90  REAL Rland_ice(niso,klon)
91  INTEGER iq
92#endif
93
94  ! FH1D
95  !     real iolat(jjm+1)
96  !real iolat(jjm+1-1/(iim*jjm))
97
98  ! Ouvrir le fichier contenant l'etat initial:
99
100  CALL open_startphy(fichnom)
101
102  ! Lecture des parametres de controle:
103
104  CALL get_var("controle", tab_cntrl)
105
106!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
108  ! Les constantes de la physiques sont lues dans la physique seulement.
109  ! Les egalites du type
110  !             tab_cntrl( 5 )=clesphy0(1)
111  ! sont remplacees par
112  !             clesphy0(1)=tab_cntrl( 5 )
113  ! On inverse aussi la logique.
114  ! On remplit les tab_cntrl avec les parametres lus dans les .def
115!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116
117  DO i = 1, length
118     tabcntr0( i ) = tab_cntrl( i )
119  ENDDO
120
121  tab_cntrl(1)=pdtphys
122  tab_cntrl(2)=radpas
123
124  ! co2_ppm : value from the previous time step
125  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
126     co2_ppm = tab_cntrl(3)
127     RCO2    = co2_ppm * 1.0e-06  * 44.011/28.97
128     ! ELSE : keep value from .def
129  END IF
130
131  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
132  co2_ppm0   = tab_cntrl(16)
133
134  solaire_etat0      = tab_cntrl(4)
135  tab_cntrl(5)=iflag_con
136  tab_cntrl(6)=nbapp_rad
137
138  if (cycle_diurne) tab_cntrl( 7) =1.
139  if (soil_model) tab_cntrl( 8) =1.
140  if (new_oliq) tab_cntrl( 9) =1.
141  if (ok_orodr) tab_cntrl(10) =1.
142  if (ok_orolf) tab_cntrl(11) =1.
143  if (ok_limitvrai) tab_cntrl(12) =1.
144
145  itau_phy = tab_cntrl(15)
146
147  clesphy0(1)=tab_cntrl( 5 )
148  clesphy0(2)=tab_cntrl( 6 )
149  clesphy0(3)=tab_cntrl( 7 )
150  clesphy0(4)=tab_cntrl( 8 )
151  clesphy0(5)=tab_cntrl( 9 )
152  clesphy0(6)=tab_cntrl( 10 )
153  clesphy0(7)=tab_cntrl( 11 )
154  clesphy0(8)=tab_cntrl( 12 )
155
156  ! set time iteration
157   CALL init_iteration(itau_phy)
158
159  ! read latitudes and make a sanity check (because already known from dyn)
160  CALL get_field("latitude",lat_startphy)
161  DO i=1,klon
162    IF (ABS(lat_startphy(i)-latitude_deg(i))>=1) THEN
163      WRITE(*,*) "phyetat0: Error! Latitude discrepancy wrt startphy file:",&
164                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
165                 " latitude_deg(i)=",latitude_deg(i)
166      ! This is presumably serious enough to abort run
167      CALL abort_physic("phyetat0","discrepancy in latitudes!",1)
168    ENDIF
169    IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.0001) THEN
170      WRITE(*,*) "phyetat0: Warning! Latitude discrepancy wrt startphy file:",&
171                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
172                 " latitude_deg(i)=",latitude_deg(i)
173    ENDIF
174  ENDDO
175
176  ! read longitudes and make a sanity check (because already known from dyn)
177  CALL get_field("longitude",lon_startphy)
178  DO i=1,klon
179    IF (ABS(lon_startphy(i)-longitude_deg(i))>=1) THEN
180      WRITE(*,*) "phyetat0: Error! Longitude discrepancy wrt startphy file:",&
181                 " i=",i," lon_startphy(i)=",lon_startphy(i),&
182                 " longitude_deg(i)=",longitude_deg(i)
183      ! This is presumably serious enough to abort run
184      CALL abort_physic("phyetat0","discrepancy in longitudes!",1)
185    ENDIF
186    IF (ABS(lon_startphy(i)-longitude_deg(i))>=0.0001) THEN
187      WRITE(*,*) "phyetat0: Warning! Longitude discrepancy wrt startphy file:",&
188                 " i=",i," lon_startphy(i)=",lon_startphy(i),&
189                 " longitude_deg(i)=",longitude_deg(i)
190    ENDIF
191  ENDDO
192
193  ! Lecture du masque terre mer
194
195  CALL get_field("masque", zmasq, found)
196  IF (.NOT. found) THEN
197     PRINT*, 'phyetat0: Le champ <masque> est absent'
198     PRINT *, 'fichier startphy non compatible avec phyetat0'
199  ENDIF
200
201  ! Lecture des fractions pour chaque sous-surface
202
203  ! initialisation des sous-surfaces
204
205  pctsrf = 0.
206
207  ! fraction de terre
208
209  CALL get_field("FTER", pctsrf(:, is_ter), found)
210  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'
211
212  ! fraction de glace de terre
213
214  CALL get_field("FLIC", pctsrf(:, is_lic), found)
215  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'
216
217  ! fraction d'ocean
218
219  CALL get_field("FOCE", pctsrf(:, is_oce), found)
220  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'
221
222  ! fraction glace de mer
223
224  CALL get_field("FSIC", pctsrf(:, is_sic), found)
225  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'
226
227  !  Verification de l'adequation entre le masque et les sous-surfaces
228
229  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
230       + pctsrf(1 : klon, is_lic)
231  DO i = 1 , klon
232     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
233        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
234             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
235             , pctsrf(i, is_lic)
236        WRITE(*, *) 'Je force la coherence zmasq=fractint'
237        zmasq(i) = fractint(i)
238     ENDIF
239  END DO
240  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
241       + pctsrf(1 : klon, is_sic)
242  DO i = 1 , klon
243     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
244        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
245             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
246             , pctsrf(i, is_sic)
247        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
248        zmasq(i) = 1. - fractint(i)
249     ENDIF
250  END DO
251
252!===================================================================
253! Lecture des temperatures du sol:
254!===================================================================
255
256  found=phyetat0_get(1,ftsol(:,1),"TS","Surface temperature",283.)
257  IF (found) THEN
258     DO nsrf=2,nbsrf
259        ftsol(:,nsrf)=ftsol(:,1)
260     ENDDO
261  ELSE
262     found=phyetat0_srf(1,ftsol,"TS","Surface temperature",283.)
263  ENDIF
264
265!===================================================================
266  ! Lecture des albedo difus et direct
267!===================================================================
268
269  DO nsrf = 1, nbsrf
270     DO isw=1, nsw
271        IF (isw.GT.99) THEN
272           PRINT*, "Trop de bandes SW"
273           call abort_physic("phyetat0", "", 1)
274        ENDIF
275        WRITE(str2, '(i2.2)') isw
276        found=phyetat0_srf(1,falb_dir(:, isw,:),"A_dir_SW"//str2//"srf","Direct Albedo",0.2)
277        found=phyetat0_srf(1,falb_dif(:, isw,:),"A_dif_SW"//str2//"srf","Direct Albedo",0.2)
278     ENDDO
279  ENDDO
280
281  found=phyetat0_srf(1,u10m,"U10M","u a 10m",0.)
282  found=phyetat0_srf(1,v10m,"V10M","v a 10m",0.)
283
284!===================================================================
285  ! Lecture des temperatures du sol profond:
286!===================================================================
287
288   DO isoil=1, nsoilmx
289        IF (isoil.GT.99) THEN
290           PRINT*, "Trop de couches "
291           call abort_physic("phyetat0", "", 1)
292        ENDIF
293        WRITE(str2,'(i2.2)') isoil
294        found=phyetat0_srf(1,tsoil(:, isoil,:),"Tsoil"//str2//"srf","Temp soil",0.)
295        IF (.NOT. found) THEN
296           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
297           PRINT*, "          Il prend donc la valeur de surface"
298           tsoil(:, isoil, :)=ftsol(:, :)
299        ENDIF
300   ENDDO
301
302!=======================================================================
303! Lecture precipitation/evaporation
304!=======================================================================
305
306  found=phyetat0_srf(1,qsurf,"QS","Near surface hmidity",0.)
307  found=phyetat0_get(1,qsol,"QSOL","Surface hmidity / bucket",0.)
308  found=phyetat0_srf(1,snow,"SNOW","Surface snow",0.)
309  found=phyetat0_srf(1,fevap,"EVAP","evaporation",0.)
310  found=phyetat0_get(1,snow_fall,"snow_f","snow fall",0.)
311  found=phyetat0_get(1,rain_fall,"rain_f","rain fall",0.)
312
313
314!=======================================================================
315! Radiation
316!=======================================================================
317
318  found=phyetat0_get(1,solsw,"solsw","net SW radiation surf",0.)
319  found=phyetat0_get(1,sollw,"sollw","net LW radiation surf",0.)
320  found=phyetat0_get(1,sollwdown,"sollwdown","down LW radiation surf",0.)
321  IF (.NOT. found) THEN
322     sollwdown = 0. ;  zts=0.
323     do nsrf=1,nbsrf
324        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
325     enddo
326     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
327  ENDIF
328
329  found=phyetat0_get(1,radsol,"RADS","Solar radiation",0.)
330  found=phyetat0_get(1,fder,"fder","Flux derivative",0.)
331
332
333  ! Lecture de la longueur de rugosite
334  found=phyetat0_srf(1,z0m,"RUG","Z0m ancien",0.001)
335  IF (found) THEN
336     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
337  ELSE
338     found=phyetat0_srf(1,z0m,"Z0m","Roughness length, momentum ",0.001)
339     found=phyetat0_srf(1,z0h,"Z0h","Roughness length, enthalpy ",0.001)
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_S",0.)
400  found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
401  found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
402  found=phyetat0_get(1,wake_fip,"WAKE_FIP","WAKE_FIP",0.)
403
404! Thermiques
405  found=phyetat0_get(1,zmax0,"ZMAX0","ZMAX0",40.)
406  found=phyetat0_get(1,f0,"F0","F0",1.e-5)
407  found=phyetat0_get(klev+1,fm_therm,"FM_THERM","Thermals mass flux",0.)
408  found=phyetat0_get(klev,entr_therm,"ENTR_THERM","Thermals Entrain.",0.)
409  found=phyetat0_get(klev,detr_therm,"DETR_THERM","Thermals Detrain.",0.)
410
411! ALE/ALP
412  found=phyetat0_get(1,ale_bl,"ALE_BL","ALE BL",0.)
413  found=phyetat0_get(1,ale_bl_trig,"ALE_BL_TRIG","ALE BL_TRIG",0.)
414  found=phyetat0_get(1,alp_bl,"ALP_BL","ALP BL",0.)
415
416!===========================================
417  ! Read and send field trs to traclmdz
418!===========================================
419
420  IF (type_trac == 'lmdz') THEN
421#ifdef ISO
422     DO it=1, nbtr         
423        iq=itr_indice(it)                                                   
424        iiq=niadv(iq)                                                           ! jyg   
425        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
426              "Surf trac"//tname(iiq),0.)
427     END DO !DO it=1, nbtr
428#else
429     DO it=1, nbtr                                                             
430!!        iiq=niadv(it+2)                                                           ! jyg
431        iiq=niadv(it+nqo)                                                           ! jyg   
432        found=phyetat0_get(1,trs(:,it),"trs_"//tname(iiq), &
433              "Surf trac"//tname(iiq),0.)
434     END DO !DO it=1, nbtr
435#endif
436     CALL traclmdz_from_restart(trs)
437
438     IF (carbon_cycle_cpl) THEN
439        ALLOCATE(co2_send(klon), stat=ierr)
440        IF (ierr /= 0) CALL abort_physic('phyetat0', 'pb allocation co2_send', 1)
441        found=phyetat0_get(1,co2_send,"co2_send","co2 send",0.)
442     END IF
443  END IF !IF (type_trac == 'lmdz') THEN
444
445#ifdef ISO
446        ! initialise les isotopes       
447        write(*,*) 'phyetat0 1069'
448         CALL phyisoetat0 (snow,run_off_lic_0, &
449     &           xtsnow,xtrun_off_lic_0, &
450     &           Rland_ice)
451#ifdef ISOVERIF
452      write(*,*) 'phyetat0 1074'
453      if (iso_eau.gt.0) then
454      call iso_verif_egalite_vect2D(  &
455     &           xtsnow,snow, &
456     &           'phyetat0 1101a',niso,klon,nbsrf)
457        do i=1,klon 
458              call iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
459     &         'phyetat0 1101b')
460         enddo
461      endif
462            write(*,*) 'phyetat0 1102'
463#endif
464#endif
465
466!===========================================
467!  ondes de gravite / relief
468!===========================================
469
470!  ondes de gravite non orographiques
471  if (ok_gwd_rando) found = &
472       phyetat0_get(klev,du_gwd_rando,"du_gwd_rando","du_gwd_rando",0.)
473  IF (.not. ok_hines .and. ok_gwd_rando) found &
474       = phyetat0_get(klev,du_gwd_front,"du_gwd_front","du_gwd_front",0.)
475
476!  prise en compte du relief sous-maille
477  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
478  found=phyetat0_get(1,zstd,"ZSTD","sub grid orography",0.)
479  found=phyetat0_get(1,zsig,"ZSIG","sub grid orography",0.)
480  found=phyetat0_get(1,zgam,"ZGAM","sub grid orography",0.)
481  found=phyetat0_get(1,zthe,"ZTHE","sub grid orography",0.)
482  found=phyetat0_get(1,zpic,"ZPIC","sub grid orography",0.)
483  found=phyetat0_get(1,zval,"ZVAL","sub grid orography",0.)
484  found=phyetat0_get(1,zmea,"ZMEA","sub grid orography",0.)
485  found=phyetat0_get(1,rugoro,"RUGSREL","sub grid orography",0.)
486
487!===========================================
488! Initialize ocean
489!===========================================
490
491  IF ( type_ocean == 'slab' ) THEN
492      CALL ocean_slab_init(dtime, pctsrf)
493      found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
494      IF (.NOT. found) THEN
495          PRINT*, "phyetat0: Le champ <tslab> est absent"
496          PRINT*, "Initialisation a tsol_oce"
497          DO i=1,nslay
498              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
499          END DO
500      END IF
501
502      ! Sea ice variables
503      found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
504      IF (version_ocean == 'sicINT') THEN
505          IF (.NOT. found) THEN
506              PRINT*, "phyetat0: Le champ <tice> est absent"
507              PRINT*, "Initialisation a tsol_sic"
508                  tice(:)=ftsol(:,is_sic)
509          END IF
510          IF (.NOT. found) THEN
511              PRINT*, "phyetat0: Le champ <seaice> est absent"
512              PRINT*, "Initialisation a 0/1m suivant fraction glace"
513              seaice(:)=0.
514              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
515                  seaice=917.
516              END WHERE
517          END IF
518      END IF !sea ice INT
519  END IF ! Slab       
520
521  ! on ferme le fichier
522  CALL close_startphy
523
524  ! Initialize module pbl_surface_mod
525
526  CALL pbl_surface_init(fder, snow, qsurf, tsoil &
527#ifdef ISO       
528     &   ,xtsnow,Rland_ice &
529#endif
530     &  )
531
532  ! Initialize module ocean_cpl_mod for the case of coupled ocean
533  IF ( type_ocean == 'couple' ) THEN
534     CALL ocean_cpl_init(dtime, longitude_deg, latitude_deg)
535  ENDIF
536
537  CALL init_iophy_new(latitude_deg, longitude_deg)
538
539  ! Initilialize module fonte_neige_mod     
540  CALL fonte_neige_init(run_off_lic_0 &
541#ifdef ISO
542    &  ,xtrun_off_lic_0 &
543#endif     
544    &   )
545
546END SUBROUTINE phyetat0
547
548!===================================================================
549FUNCTION phyetat0_get(nlev,field,name,descr,default)
550!===================================================================
551! Lecture d'un champ avec contrôle
552! Function logique dont le resultat indique si la lecture
553! s'est bien passée
554! On donne une valeur par defaut dans le cas contraire
555!===================================================================
556
557USE iostart, ONLY : get_field
558USE dimphy, only: klon
559USE print_control_mod, ONLY: lunout
560
561IMPLICIT NONE
562
563LOGICAL phyetat0_get
564
565! arguments
566INTEGER,INTENT(IN) :: nlev
567CHARACTER*(*),INTENT(IN) :: name,descr
568REAL,INTENT(IN) :: default
569REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field
570
571! Local variables
572LOGICAL found
573
574   CALL get_field(name, field, found)
575   IF (.NOT. found) THEN
576     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
577     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
578     field(:,:)=default
579   ENDIF
580   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
581   phyetat0_get=found
582
583RETURN
584END FUNCTION phyetat0_get
585
586!================================================================
587FUNCTION phyetat0_srf(nlev,field,name,descr,default)
588!===================================================================
589! Lecture d'un champ par sous-surface avec contrôle
590! Function logique dont le resultat indique si la lecture
591! s'est bien passée
592! On donne une valeur par defaut dans le cas contraire
593!===================================================================
594
595USE iostart, ONLY : get_field
596USE dimphy, only: klon
597USE indice_sol_mod, only: nbsrf
598USE print_control_mod, ONLY: lunout
599
600IMPLICIT NONE
601
602LOGICAL phyetat0_srf
603! arguments
604INTEGER,INTENT(IN) :: nlev
605CHARACTER*(*),INTENT(IN) :: name,descr
606REAL,INTENT(IN) :: default
607REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field
608
609! Local variables
610LOGICAL found,phyetat0_get
611INTEGER nsrf
612CHARACTER*2 str2
613 
614     IF (nbsrf.GT.99) THEN
615        WRITE(lunout,*) "Trop de sous-mailles"
616        call abort_physic("phyetat0", "", 1)
617     ENDIF
618
619     DO nsrf = 1, nbsrf
620        WRITE(str2, '(i2.2)') nsrf
621        found= phyetat0_get(nlev,field(:,:, nsrf), &
622        name//str2,descr//" srf:"//str2,default)
623     ENDDO
624
625     phyetat0_srf=found
626
627RETURN
628END FUNCTION phyetat0_srf
629
Note: See TracBrowser for help on using the repository browser.