source: LMDZ6/trunk/libf/phylmd/phyetat0.F90 @ 3741

Last change on this file since 3741 was 3581, checked in by oboucher, 5 years ago

Big update to the interactive carbon cycle
from Patricia's code

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