Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

limit_read_mod.F90 @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 18.0 KB

Line
1	!
2	! $Header$
3	!
4	MODULE limit_read_mod
5	!
6	! This module reads the fichier "limit.nc" containing fields for surface forcing.
7	!
8	! Module subroutines :
9	! limit_read_frac : call limit_read_tot and return the fractions
10	! limit_read_rug_alb : return rugosity and albedo, if coupled ocean call limit_read_tot first
11	! limit_read_sst : return sea ice temperature
12	! limit_read_tot : read limit.nc and store the fields in local modules variables
13	!
14	IMPLICIT NONE
15
16	REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf
17	!$OMP THREADPRIVATE(pctsrf)
18	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: rugos
19	!$OMP THREADPRIVATE(rugos)
20	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: albedo
21	!$OMP THREADPRIVATE(albedo)
22	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: sst
23	!$OMP THREADPRIVATE(sst)
24	#ifdef ISO
25	REAL, ALLOCATABLE, DIMENSION(:), SAVE, PRIVATE :: tuoce
26	!$OMP THREADPRIVATE(tuoce)
27	#endif
28	LOGICAL,SAVE :: read_continents=.FALSE.
29	!$OMP THREADPRIVATE(read_continents)
30
31	CONTAINS
32	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33	!!
34	!! Public subroutines :
35	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36
37	SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
38	!
39	! This subroutine is called from "change_srf_frac" for case of
40	! ocean=force or from ocean_slab_frac for ocean=slab.
41	! The fraction for all sub-surfaces at actual time step is returned.
42
43	USE dimphy
44	USE indice_sol_mod
45
46	! Input arguments
47	!****************************************************************************************
48	INTEGER, INTENT(IN) :: itime ! time step
49	INTEGER, INTENT(IN) :: jour ! current day
50	REAL , INTENT(IN) :: dtime ! length of time step
51
52	! Output arguments
53	!****************************************************************************************
54	REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new ! sub surface fractions
55	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
56
57	! End declaration
58	!****************************************************************************************
59
60	! 1) Read file limit.nc
61	CALL limit_read_tot(itime, dtime, jour, is_modified)
62
63	! 2) Return the fraction read in limit_read_tot
64	pctsrf_new(:,:) = pctsrf(:,:)
65
66	END SUBROUTINE limit_read_frac
67
68	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69
70	SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
71	knon, knindex, &
72	rugos_out, alb_out &
73	#ifdef ISO
74	& ,P_vegetation &
75	#endif
76	& )
77	!
78	! This subroutine is called from surf_land_bucket.
79	! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
80	! then this routine will call limit_read_tot.
81	!
82	USE dimphy
83	USE surface_data
84	#ifdef ISO
85	USE infotrac_phy, ONLY: niso
86	USE isotopes_mod, ONLY: P_veg
87	#endif
88
89	! Input arguments
90	!****************************************************************************************
91	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
92	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
93	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
94	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
95	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
96	! Output arguments
97	!****************************************************************************************
98	REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
99	REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
100
101	! Local variables
102	!****************************************************************************************
103	INTEGER :: i
104	LOGICAL :: is_modified
105
106
107	#ifdef ISO
108	real, intent(out), dimension(klon) :: P_vegetation
109	#endif
110	!****************************************************************************************
111
112	IF (type_ocean == 'couple'.OR. &
113	(type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
114	! limit.nc has not yet been read. Do it now!
115	CALL limit_read_tot(itime, dtime, jour, is_modified)
116	END IF
117
118	DO i=1,knon
119	rugos_out(i) = rugos(knindex(i))
120	alb_out(i) = albedo(knindex(i))
121	END DO
122
123	#ifdef ISO
124	!write(,) 'limit_read 120'
125	!* lecture de la fraction de l'évap sous forme de transpiration
126	!1) premier cas: on fixe P_vegetation à une valeur constante
127	! write(,) 'interfsurf 2990: P_veg=',P_veg
128	do i=1,knon
129	P_vegetation(i)=P_veg
130	enddo
131	! 2) 2e cas: on lit une carte
132	! A FAIRE
133
134	! end verif
135	#endif
136
137	END SUBROUTINE limit_read_rug_alb
138
139	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
140
141	SUBROUTINE limit_read_sst(knon, knindex, sst_out &
142	#ifdef ISO
143	& ,Roce,rlat &
144	#endif
145	)
146	!
147	! This subroutine returns the sea surface temperature already read from limit.nc.
148	! On rajoute la sortie tritium ocean sur le meme modele que la sst
149	!
150	USE dimphy, ONLY : klon
151	#ifdef ISO
152	USE infotrac_phy, ONLY: niso
153	USE isotopes_mod, ONLY: tcorr,toce,modif_sst, &
154	& deltaTtest,sstlatcrit,deltaTtestpoles,dsstlatcrit, &
155	& iso_HTO,ok_prod_nucl_tritium
156	#ifdef ISOVERIF
157	USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_positif, &
158	iso_verif_positif_nostop
159	#endif
160	#endif
161
162	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
163	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
164	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out
165	#ifdef ISO
166	REAL, DIMENSION(klon) :: tuoce_out ! sortie tritium surface ocean
167	#endif
168
169	INTEGER :: i
170	#ifdef ISO
171	real, intent(out), dimension(niso,klon) :: Roce
172	integer :: ixt
173	real, intent(in),dimension(klon) :: rlat
174	real lat_locale
175	#endif
176	!#ifdef ISOVERIF
177	! integer iso_verif_positif_nostop
178	!#endif
179
180	DO i = 1, knon
181	sst_out(i) = sst(knindex(i))
182	END DO
183
184	#ifdef ISO
185	if (iso_HTO.gt.0) then
186	if (ok_prod_nucl_tritium) then ! si on active la production nucleaire de tritium
187	DO i = 1, knon
188	tuoce_out(i)=tuoce(knindex(i))
189	END DO
190	endif
191	endif
192	#endif
193
194	#ifdef ISO
195	if (modif_sst.ge.1) then
196	do i = 1, knon
197	lat_locale=rlat(knindex(i))
198	! test: modification uniforme de la sst
199	if (modif_sst.eq.1) then
200	sst_out(i)= sst_out(i)+deltaTtest
201	elseif (modif_sst.eq.2) then !if (modif_sst.eq.1) then
202	! pattern parabolique en dehors des tropiques (sstlatcrit)
203	if (abs(lat_locale).gt.sstlatcrit) then
204	sst_out(i)= sst_out(i)+deltaTtestpoles &
205	& (lat_locale2-sstlatcrit*2) &
206	& /(90.02-sstlatcrit2)
207	endif !if (abs(lat_locale).gt.abs(sstlatcrit)) then
208
209	else if (modif_sst.eq.3) then
210
211	if (abs(lat_locale).gt.abs(sstlatcrit)) then
212	if (abs(lat_locale).gt.sstlatcrit+dsstlatcrit) then
213	sst_out(i)= sst_out(i)+deltaTtestpoles
214	else
215	sst_out(i)= sst_out(i)+deltaTtestpoles &
216	& *(abs(lat_locale)-sstlatcrit)/dsstlatcrit
217	endif
218	endif
219	endif !if (modif_sst.eq.1) then
220	enddo !do i = 1, knon
221	endif !if (modif_sst.ge.1) then
222	#endif
223	#ifdef ISOVERIF
224	do i=1,knon,20
225	call iso_verif_positif(sst_out(i)-100.0,'limit_read 4323')
226	enddo
227	#endif
228
229
230	#ifdef ISO
231	!* lecture de Roce
232	! 1) première possibilité: valeur fixe à SMOW
233	DO i = 1, knon
234	do ixt=1,niso
235	Roce(ixt,i)=tcorr(ixt)*toce(ixt)
236	enddo !do ixt=1,niso
237	enddo !DO i = 1, knon
238	! 2) deuxième possibilité: lecture de la carte
239	! A FAIRE
240
241	! lecture pour le tritium
242	if ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) then
243	! lecture de la carte tritium ocean surface
244	Roce(iso_HTO,i)=tcorr(iso_HTO)tuoce_out(i)1.E-18*2.
245	endif
246	#endif
247
248	#ifdef ISOVERIF
249	do i=1,knon
250	if (iso_verif_positif_nostop(370.0-sst_out(i), &
251	'limit_read 368').eq.1) then
252	write(,) 'i,knindex,sst_out=',i,knindex,sst_out(i)
253	stop
254	endif
255	enddo
256	#endif
257
258	END SUBROUTINE limit_read_sst
259
260	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261	!!
262	!! Private subroutine :
263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264
265	SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
266	!
267	! Read everything needed from limit.nc
268	!
269	! 0) Initialize
270	! 1) Open the file limit.nc, if it is time
271	! 2) Read fraction, if not type_ocean=couple
272	! 3) Read sea surface temperature, if not type_ocean=couple
273	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
274	! 5) Close file and distribuate variables to all processus
275
276	USE dimphy
277	USE mod_grid_phy_lmdz
278	USE mod_phys_lmdz_para
279	USE surface_data, ONLY : type_ocean, ok_veget
280	USE netcdf
281	USE indice_sol_mod
282	#ifdef ISO
283	!USE infotrac_phy, ONLY: use_iso
284	USE isotopes_mod, ONLY : iso_HTO,ok_prod_nucl_tritium
285	#ifdef ISOVERIF
286	USE isotopes_verif_mod, ONLY : iso_verif_positif_nostop
287	#endif
288	#endif
289
290	IMPLICIT NONE
291
292	! In- and ouput arguments
293	!****************************************************************************************
294	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
295	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
296	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
297
298	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
299
300	! Locals variables with attribute SAVE
301	!****************************************************************************************
302	! frequence de lecture des conditions limites (en pas de physique)
303	INTEGER,SAVE :: lmt_pas
304	!$OMP THREADPRIVATE(lmt_pas)
305	LOGICAL, SAVE :: first_call=.TRUE.
306	!$OMP THREADPRIVATE(first_call)
307	INTEGER, SAVE :: jour_lu = -1
308	!$OMP THREADPRIVATE(jour_lu)
309	! Locals variables
310	!****************************************************************************************
311	INTEGER :: nid, nvarid
312	INTEGER :: ii, ierr
313	INTEGER, DIMENSION(2) :: start, epais
314	REAL, DIMENSION(klon_glo,nbsrf) :: pct_glo ! fraction at global grid
315	REAL, DIMENSION(klon_glo) :: sst_glo ! sea-surface temperature at global grid
316	REAL, DIMENSION(klon_glo) :: rug_glo ! rugosity at global grid
317	REAL, DIMENSION(klon_glo) :: alb_glo ! albedo at global grid
318	CHARACTER(len=20) :: modname='limit_read_mod'
319	#ifdef ISO
320	REAL, DIMENSION(klon_glo) :: tuoce_glo ! sea-surface tritium et global grid
321	#endif
322
323	! End declaration
324	!****************************************************************************************
325
326	!****************************************************************************************
327	! 0) Initialization
328	!
329	!****************************************************************************************
330	IF (first_call) THEN
331	! calculate number of time steps for one day
332	lmt_pas = NINT(86400./dtime * 1.0)
333
334	! Allocate module save variables
335	IF ( type_ocean /= 'couple' ) THEN
336	ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
337	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
338	END IF
339
340	IF ( .NOT. ok_veget ) THEN
341	ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
342	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
343	END IF
344
345	first_call=.FALSE.
346	ENDIF
347
348	!****************************************************************************************
349	! 1) Open the file limit.nc if it is the right moment to read, once a day.
350	! The file is read only by the master thread of the master mpi process(is_mpi_root)
351	!
352	!****************************************************************************************
353
354	is_modified = .FALSE.
355	IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN ! time to read
356	jour_lu = jour
357	is_modified = .TRUE.
358	!$OMP MASTER ! Only master thread
359	IF (is_mpi_root) THEN ! Only master processus
360
361	ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
362	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
363	'Pb d''ouverture du fichier de conditions aux limites',1)
364
365	! La tranche de donnees a lire:
366	start(1) = 1
367	start(2) = jour
368	epais(1) = klon_glo
369	epais(2) = 1
370
371
372	!****************************************************************************************
373	! 2) Read fraction if not type_ocean=couple
374	!
375	!****************************************************************************************
376
377	IF ( type_ocean /= 'couple') THEN
378	!
379	! Ocean fraction
380	ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
381	IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
382
383	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
384	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
385	!
386	! Sea-ice fraction
387	ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
388	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
389
390	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
391	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
392
393
394	! Read land and continentals fraction only if asked for
395	IF (read_continents .OR. itime == 1) THEN
396	!
397	! Land fraction
398	ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
399	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
400
401	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
402	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
403	!
404	! Continentale ice fraction
405	ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
406	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
407
408	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
409	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
410	END IF
411
412	END IF ! type_ocean /= couple
413
414	!****************************************************************************************
415	! 3) Read sea-surface temperature, if not coupled ocean
416	!
417	!****************************************************************************************
418	IF ( type_ocean /= 'couple') THEN
419
420	ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
421	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
422
423	ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
424	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
425
426	#ifdef ISO
427	IF ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) THEN
428	ierr = NF90_INQ_VARID(nid, 'TUOCE', nvarid)
429	IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Le champ <TUOCE> est absent',1)
430
431	ierr = NF90_GET_VAR(nid,nvarid,tuoce_glo,start,epais)
432	IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Lecture echouee pour <TUOCE>',1)
433	END IF
434	#ifdef ISOVERIF
435	do ii=1,klon_glo
436	if (iso_verif_positif_nostop(370.0-sst_glo(ii), &
437	'limit_read 384').eq.1) then
438	write(,) 'ii,sst_glo=',ii,sst_glo(ii)
439	write(,) 'jour,start,epais=',jour,start,epais
440	stop
441	endif
442	enddo
443	#endif
444	#endif
445
446	END IF
447
448	!****************************************************************************************
449	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
450	!
451	!****************************************************************************************
452
453	IF (.NOT. ok_veget) THEN
454	!
455	! Read albedo
456	ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
457	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
458
459	ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
460	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
461	!
462	! Read rugosity
463	ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
464	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
465
466	ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
467	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
468
469	END IF
470
471	!****************************************************************************************
472	! 5) Close file and distribuate variables to all processus
473	!
474	!****************************************************************************************
475	ierr = NF90_CLOSE(nid)
476	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
477	ENDIF ! is_mpi_root
478
479	!$OMP END MASTER
480	!$OMP BARRIER
481
482	IF ( type_ocean /= 'couple') THEN
483	CALL Scatter(sst_glo,sst)
484	CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
485	CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
486	IF (read_continents .OR. itime == 1) THEN
487	CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
488	CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
489	END IF
490	#ifdef ISO
491	IF ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) THEN
492	CALL Scatter(tuoce_glo,tuoce)
493	END IF
494	#endif
495	END IF
496
497	IF (.NOT. ok_veget) THEN
498	CALL Scatter(alb_glo, albedo)
499	CALL Scatter(rug_glo, rugos)
500	END IF
501
502	ENDIF ! time to read
503
504	END SUBROUTINE limit_read_tot
505
506
507	END MODULE limit_read_mod

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

Download in other formats:

Original Format