Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

icolmdz_etat0.f90 @ 5799

Last change on this file since 5799 was 5590, checked in by yann meurdesoif, 8 months ago
Add phylmdiso directory for testing ICOLMDZISO configuration. YM
Property svn:eol-style set to `native` Property svn:executable set to ``*
File size: 16.2 KB

Line
1	MODULE icolmdz_etat0
2
3
4	CONTAINS
5
6	SUBROUTINE init_etat0
7	USE xios_mod
8	USE omp_para
9	USE getin_mod, ONLY: getin
10	IMPLICIT NONE
11	CHARACTER(LEN=256) :: etat0_lmdz
12
13	! @getin_doc etat0_database_type realm=initial type=String default=`legacy`
14	! When `etat0=="database"`, selects type of external file(s) from which initial state is read. Valid values :
15	! * `legacy` : reads from XIOS group `read_files_legacy`, see `xml/context_dynamico.xml`
16	! * `ERA5_forcing` : reads from XIOS group `read_files_ERA5_forcing`, see `xml/context_dynamico.xml`
17
18	etat0_lmdz = "legacy"
19	CALL getin("etat0_lmdz", etat0_lmdz)
20
21	IF (is_omp_master) THEN
22	CALL xios_set_fieldgroup_attr("read_fields", enabled=.TRUE.)
23	IF (TRIM(etat0_lmdz) == "legacy") THEN
24	CALL xios_set_filegroup_attr("read_files_legacy", enabled=.TRUE.)
25	CALL xios_set_field_attr("relief_db", field_ref="relief_db_legacy")
26	CALL xios_set_field_attr("ps_db", field_ref="ps_db_legacy")
27	CALL xios_set_field_attr("z_db", field_ref="z_db_legacy")
28	CALL xios_set_field_attr("ts_db", field_ref="ts_db_legacy")
29	CALL xios_set_field_attr("u_db", field_ref="u_db_legacy")
30	CALL xios_set_field_attr("v_db", field_ref="v_db_legacy")
31	CALL xios_set_field_attr("temp_db", field_ref="temp_db_legacy")
32	CALL xios_set_field_attr("q_db", field_ref="q_db_legacy")
33	ELSE IF (TRIM(etat0_lmdz) == "ERA5_forcing") THEN
34	CALL xios_set_filegroup_attr("read_files_ERA5_forcing", enabled=.TRUE.)
35	CALL xios_set_field_attr("relief_db", field_ref="relief_db_forcing")
36	CALL xios_set_field_attr("ps_db", field_ref="ps_db_forcing")
37	CALL xios_set_field_attr("z_db", field_ref="z_db_forcing")
38	CALL xios_set_field_attr("ts_db", field_ref="ts_db_forcing")
39	CALL xios_set_field_attr("u_db", field_ref="u_db_forcing")
40	CALL xios_set_field_attr("v_db", field_ref="v_db_forcing")
41	CALL xios_set_field_attr("temp_db", field_ref="temp_db_forcing")
42	CALL xios_set_field_attr("q_db", field_ref="q_db_forcing")
43	ELSE
44	PRINT*,"Bad selector for variable <etat0_lmdz> ",TRIM(etat0_lmdz), " option are <legacy> (default), <ERA5_forcing>"
45	STOP
46	END IF
47	END IF
48	END SUBROUTINE init_etat0
49
50
51	SUBROUTINE etat0(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
52	USE icosa
53	USE restart_mod
54	USE wind_from_lonlat_mod
55	USE write_field_mod
56	USE time_mod
57	USE transfert_mod
58	USE xios_mod
59	USE write_field_mod
60	USE vertical_remap_mod
61	USE theta2theta_rhodz_mod
62	USE qsat_mod
63	USE pression_mod
64	USE omp_para
65	USE tracer_icosa_mod
66	IMPLICIT NONE
67	TYPE(t_field), POINTER :: f_ps(:)
68	TYPE(t_field), POINTER :: f_phis(:)
69	TYPE(t_field), POINTER :: f_theta_rhodz(:)
70	TYPE(t_field), POINTER :: f_u(:)
71	TYPE(t_field), POINTER :: f_q(:)
72
73	TYPE(t_field), POINTER, SAVE :: f_ulon_reg(:)
74	TYPE(t_field), POINTER, SAVE :: f_ulat_reg(:)
75	TYPE(t_field), POINTER, SAVE :: f_temp_reg(:)
76	TYPE(t_field), POINTER, SAVE :: f_q_reg(:)
77
78	TYPE(t_field), POINTER, SAVE :: f_ts(:)
79	TYPE(t_field), POINTER, SAVE :: f_z(:)
80	TYPE(t_field), POINTER, SAVE :: f_ulon(:)
81	TYPE(t_field), POINTER, SAVE :: f_ulat(:)
82	TYPE(t_field), POINTER, SAVE :: f_temp(:)
83	TYPE(t_field), POINTER, SAVE :: f_q1(:)
84	TYPE(t_field), POINTER, SAVE :: f_qsat(:)
85	TYPE(t_field), POINTER, SAVE :: f_p(:)
86	INTEGER :: nb_level
87	REAL, ALLOCATABLE:: levels(:)
88	INTEGER :: ind
89	INTEGER :: iq, iq_H2Og, iq_H2Ol, iq_H2Os
90
91	CALL xios_read_field("relief_db", f_phis)
92
93	CALL writeField("relief_out", f_phis, once=.TRUE.)
94
95	IF (is_omp_level_master) THEN
96	DO ind = 1, ndomain
97	IF (.NOT. assigned_domain(ind)) CYCLE
98	f_phis(ind)%rval2d(:) = f_phis(ind)%rval2d(:)*g
99	END DO
100	END IF
101	!$OMP BARRIER
102
103	IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", n_glo=nb_level)
104	CALL bcast_omp(nb_level)
105	ALLOCATE (levels(nb_level))
106
107	IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", value=levels)
108	CALL bcast_omp(levels)
109
110	levels = levels*100 ! hectoPascal -> Pascal
111
112	CALL allocate_field(f_ts, field_t, type_real, name="ts")
113	CALL allocate_field(f_z, field_t, type_real, name="z")
114	CALL allocate_field(f_ulon_reg, field_t, type_real, nb_level)
115	CALL allocate_field(f_ulat_reg, field_t, type_real, nb_level)
116	CALL allocate_field(f_temp_reg, field_t, type_real, nb_level)
117	CALL allocate_field(f_q_reg, field_t, type_real, nb_level)
118
119	CALL allocate_field(f_q1, field_t, type_real, llm)
120	CALL allocate_field(f_qsat, field_t, type_real, llm)
121	CALL allocate_field(f_p, field_t, type_real, llm + 1)
122	CALL allocate_field(f_temp, field_t, type_real, llm)
123	CALL allocate_field(f_ulon, field_t, type_real, llm)
124	CALL allocate_field(f_ulat, field_t, type_real, llm)
125
126	CALL xios_read_field("z_db", f_z)
127	CALL xios_read_field("ps_db", f_ps)
128	CALL xios_read_field("ts_db", f_ts)
129	CALL writeField("ps_out", f_ps)
130
131	!$OMP BARRIER
132
133	! CALL writeField("phis_out",f_phis,once=.TRUE.)
134	! CALL writeField("ts_out",f_ts,once=.TRUE.)
135
136	! make correction to ps due to relief at higher resolution
137	! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY
138	IF (is_omp_level_master) THEN
139	DO ind = 1, ndomain
140	IF (.NOT. assigned_domain(ind)) CYCLE
141	f_ps(ind)%rval2d(:) = f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:) - f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:))
142	END DO
143	END IF
144	!$OMP BARRIER
145	CALL transfert_request(f_ps, req_i0)
146	CALL writeField("ps_out", f_ps)
147
148	CALL xios_read_field("temp_db", f_temp_reg)
149	CALL vertical_remap(levels, f_temp_reg, f_ps, f_temp)
150	CALL transfert_request(f_temp, req_i0)
151	CALL temperature2theta_rhodz(f_ps, f_temp, f_theta_rhodz)
152
153	CALL xios_read_field("u_db", f_ulon_reg)
154	CALL vertical_remap(levels, f_ulon_reg, f_ps, f_ulon)
155	CALL xios_read_field("v_db", f_ulat_reg)
156	CALL vertical_remap(levels, f_ulat_reg, f_ps, f_ulat)
157	CALL transfert_request(f_ulat, req_i0)
158	CALL transfert_request(f_ulon, req_i0)
159	CALL ulonlat2un(f_ulon, f_ulat, f_u)
160
161	CALL xios_read_field("q_db", f_q_reg)
162	CALL vertical_remap(levels, f_q_reg, f_ps, f_q1)
163
164	CALL pression(f_ps, f_p)
165	! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2
166	CALL qsat(f_temp, f_p, f_qsat)
167	CALL transfert_request(f_qsat, req_i0)
168
169	DO iq = 1, nqtot
170	IF (TRIM(tracers(iq)%name) == 'H2O_g') iq_H2Og=iq
171	IF (TRIM(tracers(iq)%name) == 'H2O_l') iq_H2Ol=iq
172	IF (TRIM(tracers(iq)%name) == 'H2O_s') iq_H2Os=iq
173	ENDDO
174
175	DO ind = 1, ndomain
176	IF (.NOT. assigned_domain(ind)) CYCLE
177
178	IF (tracers(iq_H2Ol)%has_default_init_value) THEN
179	! f_q(ind)%rval4d(:, :, iq_H2Ol) = tracers(iq_H2Ol)%default_init_value
180	f_q(ind)%rval4d(:, :, iq_H2Ol) = tracers(iq_H2Ol)%default_init_value
181	ELSE
182	f_q(ind)%rval4d(:, :, iq_H2Ol) = 1e-6
183	END IF
184
185	f_q(ind)%rval4d(:, :, iq_H2Og) = f_q1(ind)%rval3d(:, :)f_qsat(ind)%rval3d(:, :)0.01
186	WHERE(f_q(ind)%rval4d(:, :, iq_H2Og)<1e-20) f_q(ind)%rval4d(:, :, iq_H2Og) = 1e-20
187	f_q(ind)%rval4d(:, :, iq_H2Os) = 1e-10
188	END DO
189
190	tracers(iq_H2Og)%already_initialized = .TRUE.
191	tracers(iq_H2Ol)%already_initialized = .TRUE.
192	tracers(iq_H2Os)%already_initialized = .TRUE.
193
194
195
196	PRINT*,'CALLING ISOTOPE INITIALIZATION'
197	CALL etat0_iso(f_q)
198
199	! CALL switch_iso(f_q)
200
201	CALL writeField("tempdb_out", f_temp_reg)
202	CALL writeField("temp_out", f_temp)
203
204	CALL deallocate_field(f_ts)
205	CALL deallocate_field(f_z)
206	CALL deallocate_field(f_ulon_reg)
207	CALL deallocate_field(f_ulat_reg)
208	CALL deallocate_field(f_temp_reg)
209	CALL deallocate_field(f_q_reg)
210
211	CALL deallocate_field(f_q1)
212	CALL deallocate_field(f_qsat)
213	CALL deallocate_field(f_p)
214	CALL deallocate_field(f_temp)
215	CALL deallocate_field(f_ulon)
216	CALL deallocate_field(f_ulat)
217
218	END SUBROUTINE etat0
219
220	SUBROUTINE switch_iso(f_q)
221	USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ
222	USE isotopes_routines_mod, ONLY : fractcalk_liq, calcul_kcin ! ==> LMDZ
223	USE infotrac_phy, ONLY: isoName, niso ! ==> LMDZ
224	USE icosa ! ==> DYNAMICO
225	USE tracer_icosa_mod
226	USE field_mod
227	USE omp_para, ONLY: is_master
228	IMPLICIT NONE
229	TYPE(t_field), POINTER :: f_q(:)
230	INTEGER :: ind, iq, iq_H2Og, iq_H216Og, iq_H218Og, iq_H2Ol
231
232	DO iq = 1, nqtot
233	IF (TRIM(tracers(iq)%name) == 'H2O_l') iq_H2Ol=iq
234	IF (TRIM(tracers(iq)%name) == 'H2O_g') iq_H2Og=iq
235	IF (TRIM(tracers(iq)%name) == 'H216O_g') iq_H216Og=iq
236	IF (TRIM(tracers(iq)%name) == 'H218O_g') iq_H218Og=iq
237	ENDDO
238
239	DO ind = 1, ndomain
240	IF (.NOT. assigned_domain(ind)) CYCLE
241	f_q(ind)%rval4d(:, :, iq_H2Ol) = f_q(ind)%rval4d(:, :, iq_H2Og)
242	f_q(ind)%rval4d(:, :, iq_H218Og) = f_q(ind)%rval4d(:, :, iq_H2Og)
243	f_q(ind)%rval4d(:, :, iq_H2Og) = 1.
244	ENDDO
245
246	END SUBROUTINE switch_iso
247
248	SUBROUTINE etat0_iso(f_q)
249	USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ
250	USE isotopes_routines_mod, ONLY : fractcalk_liq, calcul_kcin ! ==> LMDZ
251	USE infotrac_phy, ONLY: isoName, niso ! ==> LMDZ
252	USE icosa ! ==> DYNAMICO
253	USE tracer_icosa_mod
254	USE field_mod
255	USE omp_para, ONLY: is_master
256	USE compute_transport_mod, ONLY : q_minimum
257	IMPLICIT NONE
258	TYPE(t_field), POINTER :: f_q(:)
259
260	INTEGER :: iq, iqchild, iq_wiso
261	CHARACTER(LEN=256) :: wiso_name
262
263	! idealised conditions for Rayleigh distillation
264	! R = R0_evap(wiso) * (q / q0_evap)**(alpha_ideal(wiso) - 1)
265	REAL,PARAMETER :: q0_evap = 20e-3 ! initial condition of water vapour (kg/kg) : for Rayleigh distillation
266	REAL,PARAMETER :: rh0_evap = 0.7 ! initial condition of relative humidity (1) : for R0_evap
267	REAL :: ws0_evap = 2.0 ! initial condition of wind speed (m/s) : for R0_evap
268	REAL :: kcin(niso) ! kinetic coefficient at evaporation = f(ws0_evap) : for R0_evap
269	REAL,PARAMETER :: temp_ideal = 283. ! mean atmospheric temperature used for fractionation (K) : for Rayleigh distillation
270	! isotope-dependent
271	REAL :: R0_evap ! mean value of initial isotopic ratio of evaporation over ocean : for Rayleigh distillation
272	REAL :: alpha_ideal ! mean value of fractionation coefficient : for Rayleigh distillation
273	REAL :: Rocean ! mean value of initial isotopic ratio of ocean : for R0_evap (m/s)
274	REAL :: Rdefault_val ! default value of initial isotopic ratio
275	REAL :: kcin_ocean ! kinetic coefficient at evaporation = f(ws0_evap) : for R0_evap
276
277	INTEGER :: ind
278
279	PRINT*, "ISOTOPE INITIALIZED ? ", using_iso
280	IF (.NOT. using_iso) RETURN
281
282	! overwrite initialisation_iso => start from scratch
283	initialisation_iso = 1
284
285	CALL calcul_kcin(ws0_evap, kcin)
286
287	DO ind = 1, ndomain
288	IF (.NOT. assigned_domain(ind)) CYCLE
289	! CAa : loop on the tracers: look for water components
290	DO iq = 1, nqtot
291	! looking for vapour = 'H2O_g'
292	! "from readTracFiles_mod: (Only convertable names are water descendants names H2O_<phase>, <isotope>_<phase>, <isotope>_<phase>_<tag>, with:)"
293	IF (TRIM(tracers(iq)%name) == 'H2O_g') THEN
294	! Loop on children of water vapour
295	! ================================
296	! idealised Rayleigh distillation
297
298	DO iqchild = 1, tracers(iq)%nb_children
299	iq_wiso = tracers(iq)%children(iqchild)
300	wiso_name = tracers(iq_wiso)%name
301	SELECT CASE(TRIM(wiso_name))
302	CASE('H216O_g');
303	CALL fractcalk_liq(iso_eau, temp_ideal, alpha_ideal)
304	Rocean = toce(iso_eau)
305	kcin_ocean= kcin(iso_eau)
306	CASE('H217O_g')
307	CALL fractcalk_liq(iso_O17, temp_ideal, alpha_ideal)
308	Rocean = toce(iso_O17)
309	kcin_ocean= kcin(iso_O17)
310	CASE('H218O_g')
311	CALL fractcalk_liq(iso_O18, temp_ideal, alpha_ideal)
312	Rocean = toce(iso_O18)
313	kcin_ocean= kcin(iso_O18)
314	CASE('HDO_g')
315	CALL fractcalk_liq(iso_HDO, temp_ideal, alpha_ideal)
316	Rocean = toce(iso_HDO)
317	kcin_ocean= kcin(iso_O18)
318	CASE('HTO_g')
319	CALL fractcalk_liq(iso_HTO, temp_ideal, alpha_ideal)
320	Rocean = toce(iso_HTO)
321	kcin_ocean= kcin(iso_O18)
322	CASE DEFAULT
323	PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP'
324	STOP
325	END SELECT
326
327	! R0_evap = R evaporation from Merlivat
328	R0_evap = Rocean / alpha_ideal * (1.0 - kcin_ocean) / (1.0 - kcin_ocean * rh0_evap)
329	! R_iso = R from idealised Rayleigh distillation
330	! R_iso = R0_evap(wiso) * (q / q0_evap)**(alpha_ideal(wiso) - 1)
331	f_q(ind)%rval4d(:, :, iq_wiso) = R0_evap * (min(q0_evap, f_q(ind)%rval4d(:, :, iq)) / q0_evap)**(alpha_ideal - 1)
332	IF (is_master) PRINT*,TRIM(wiso_name)," ",Rocean, kcin_ocean, R0_evap, alpha_ideal, MAXVAL(f_q(ind)%rval4d(:, :, iq)),MINVAL(f_q(ind)%rval4d(:, :, iq)), MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso))
333	tracers(iq_wiso)%already_initialized = .TRUE.
334	END DO
335
336	ELSE IF (TRIM(tracers(iq)%name) == 'H2O_l') THEN
337	! TODO : update, it's wrong
338
339	DO iqchild = 1, tracers(iq)%nb_children
340	iq_wiso = tracers(iq)%children(iqchild)
341	wiso_name = tracers(iq_wiso)%name
342	SELECT CASE(wiso_name)
343	CASE('H216O_l')
344	Rdefault_val = Rdefault(iso_eau)
345	CASE('H217O_l')
346	Rdefault_val = Rdefault(iso_O17)
347	CASE('H218O_l')
348	Rdefault_val = Rdefault(iso_O18)
349	CASE('HDO_l')
350	Rdefault_val = Rdefault(iso_HDO)
351	CASE('HTO_l')
352	Rdefault_val = Rdefault(iso_HTO)
353	CASE DEFAULT
354	PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP'
355	STOP
356	END SELECT
357	f_q(ind)%rval4d(:, :, iq_wiso) = Rdefault_val
358	IF (is_master) PRINT*,TRIM(wiso_name)," ",Rdefault_val, MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso))
359	tracers(iq_wiso)%already_initialized = .TRUE.
360	END DO
361
362	ELSE IF (TRIM(tracers(iq)%name) == 'H2O_s') THEN
363
364	DO iqchild = 1, tracers(iq)%nb_children
365	iq_wiso = tracers(iq)%children(iqchild)
366	wiso_name = tracers(iq_wiso)%name
367	SELECT CASE(wiso_name)
368	CASE('H216O_s')
369	Rdefault_val = Rdefault(iso_eau)
370	CASE('H217O_s')
371	Rdefault_val = Rdefault(iso_O17)
372	CASE('H218O_s')
373	Rdefault_val = Rdefault(iso_O18)
374	CASE('HDO_s')
375	Rdefault_val = Rdefault(iso_HDO)
376	CASE('HTO_s')
377	Rdefault_val = Rdefault(iso_HTO)
378	CASE DEFAULT
379	PRINT*, 'icosa_lmdz isotope init"' // TRIM(wiso_name) // '" does not exists, STOP'
380	STOP
381	END SELECT
382	f_q(ind)%rval4d(:, :, iq_wiso) = Rdefault_val
383	IF (is_master) PRINT*,TRIM(wiso_name)," ",Rdefault_val, MAXVAL(f_q(ind)%rval4d(:, :, iq_wiso)),MINVAL(f_q(ind)%rval4d(:, :, iq_wiso))
384	tracers(iq_wiso)%already_initialized = .TRUE.
385	END DO
386	END IF
387
388	ENDDO
389
390	ENDDO
391
392
393	END SUBROUTINE etat0_iso
394
395
396	END MODULE icolmdz_etat0

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

Download in other formats:

Original Format