Context Navigation

phyetat0.F90 @ 4298

Last change on this file since 4298 was 4298, checked in by pcadule, 21 months ago
modifications to ensure restartability of the model when CO2 tracer is passed to radiative code, and to add diagnostics variables
File size: 25.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format