source: LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90 @ 4358

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