source: LMDZ6/trunk/libf/phylmdiso/phyetat0.F90 @ 4071

Last change on this file since 4071 was 4071, checked in by dcugnet, 2 years ago
  • Fix for unadvected tracers (iadv==0)
  • The key %isH2Ofamily, from the derived type "trac_type", is replaced with the more general

key %isInPhysics, which is TRUE for tracers both in "qx" and "tr_seri".

Currently, FALSE for tracers descending on H2O (isotopes and tagging tracers included). Could be set to FALSE
for interactive CO2 (type_trac=='inco') or ice supersaturated cloud content (tranfered to "rneb_seri")

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