Context Navigation

limit_read_mod.f90 @ 5644

Last change on this file since 5644 was 5549, checked in by aborella, 5 months ago
Fix for reading aviation data with multiple OMPs
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.8 KB

Line
1	!
2	! $Id: limit_read_mod.f90 5549 2025-02-18 16:49:06Z aborella $
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	LOGICAL,SAVE :: read_continents=.FALSE.
25	!$OMP THREADPRIVATE(read_continents)
26
27	CONTAINS
28	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29	!!
30	!! Public subroutines :
31	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32
33
34	SUBROUTINE init_limit_read(first_day)
35	USE mod_grid_phy_lmdz
36	USE surface_data
37	USE mod_phys_lmdz_para
38	USE lmdz_xios
39
40	IMPLICIT NONE
41	INTEGER, INTENT(IN) :: first_day
42
43	IF ( type_ocean /= 'couple') THEN
44	IF (grid_type==unstructured) THEN
45	IF (is_omp_master) CALL xios_set_file_attr("limit_read",enabled=.TRUE.,record_offset=first_day)
46	ENDIF
47	ENDIF
48
49	END SUBROUTINE init_limit_read
50
51	SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
52	!
53	! This subroutine is called from "change_srf_frac" for case of
54	! ocean=force or from ocean_slab_frac for ocean=slab.
55	! The fraction for all sub-surfaces at actual time step is returned.
56
57	USE dimphy
58	USE indice_sol_mod
59
60	! Input arguments
61	!****************************************************************************************
62	INTEGER, INTENT(IN) :: itime ! time step
63	INTEGER, INTENT(IN) :: jour ! current day
64	REAL , INTENT(IN) :: dtime ! length of time step
65
66	! Output arguments
67	!****************************************************************************************
68	REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new ! sub surface fractions
69	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
70
71	! End declaration
72	!****************************************************************************************
73
74	! 1) Read file limit.nc
75	CALL limit_read_tot(itime, dtime, jour, is_modified)
76
77	! 2) Return the fraction read in limit_read_tot
78	pctsrf_new(:,:) = pctsrf(:,:)
79
80	END SUBROUTINE limit_read_frac
81
82	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
83
84	SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
85	knon, knindex, &
86	rugos_out, alb_out)
87	!
88	! This subroutine is called from surf_land_bucket.
89	! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
90	! then this routine will call limit_read_tot.
91	!
92	USE dimphy
93	USE surface_data
94
95	! Input arguments
96	!****************************************************************************************
97	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
98	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
99	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
100	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
101	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
102	! Output arguments
103	!****************************************************************************************
104	REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
105	REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
106
107	! Local variables
108	!****************************************************************************************
109	INTEGER :: i
110	LOGICAL :: is_modified
111	!****************************************************************************************
112
113	IF (type_ocean == 'couple'.OR. &
114	(type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
115	! limit.nc has not yet been read. Do it now!
116	CALL limit_read_tot(itime, dtime, jour, is_modified)
117	END IF
118
119	DO i=1,knon
120	rugos_out(i) = rugos(knindex(i))
121	alb_out(i) = albedo(knindex(i))
122	END DO
123
124	END SUBROUTINE limit_read_rug_alb
125
126	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
127
128	SUBROUTINE limit_read_sst(knon, knindex, sst_out)
129	!
130	! This subroutine returns the sea surface temperature already read from limit.nc.
131	!
132	USE dimphy, ONLY : klon
133
134	INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid
135	INTEGER, DIMENSION(klon), INTENT(IN) :: knindex ! grid point number for compressed grid
136	REAL, DIMENSION(klon), INTENT(OUT) :: sst_out
137
138	INTEGER :: i
139
140	DO i = 1, knon
141	sst_out(i) = sst(knindex(i))
142	END DO
143
144	END SUBROUTINE limit_read_sst
145
146	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
147	!!
148	!! Private subroutine :
149	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150
151	SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
152	!
153	! Read everything needed from limit.nc
154	!
155	! 0) Initialize
156	! 1) Open the file limit.nc, if it is time
157	! 2) Read fraction, if not type_ocean=couple
158	! 3) Read sea surface temperature, if not type_ocean=couple
159	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
160	! 5) Close file and distribuate variables to all processus
161
162	USE dimphy
163	USE mod_grid_phy_lmdz
164	USE mod_phys_lmdz_para
165	USE surface_data, ONLY : type_ocean, ok_veget
166	USE netcdf
167	USE indice_sol_mod
168	USE phys_cal_mod, ONLY : calend, year_len
169	USE print_control_mod, ONLY: lunout, prt_level
170	USE lmdz_xios, ONLY: xios_recv_field, using_xios
171
172	IMPLICIT NONE
173
174	! In- and ouput arguments
175	!****************************************************************************************
176	INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
177	INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
178	REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
179
180	LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
181
182	! Locals variables with attribute SAVE
183	!****************************************************************************************
184	! frequence de lecture des conditions limites (en pas de physique)
185	INTEGER,SAVE :: lmt_pas
186	!$OMP THREADPRIVATE(lmt_pas)
187	LOGICAL, SAVE :: first_call=.TRUE.
188	!$OMP THREADPRIVATE(first_call)
189	INTEGER, SAVE :: jour_lu = -1
190	!$OMP THREADPRIVATE(jour_lu)
191	! Locals variables
192	!****************************************************************************************
193	INTEGER :: nid, nvarid, ndimid, nn
194	INTEGER :: ii, ierr
195	INTEGER, DIMENSION(2) :: start, epais
196	REAL, DIMENSION(klon_glo,nbsrf) :: pct_glo ! fraction at global grid
197	REAL, DIMENSION(klon_glo) :: sst_glo ! sea-surface temperature at global grid
198	REAL, DIMENSION(klon_glo) :: rug_glo ! rugosity at global grid
199	REAL, DIMENSION(klon_glo) :: alb_glo ! albedo at global grid
200
201	REAL, DIMENSION(klon_mpi,nbsrf) :: pct_mpi ! fraction at global grid
202	REAL, DIMENSION(klon_mpi) :: sst_mpi ! sea-surface temperature at global grid
203	REAL, DIMENSION(klon_mpi) :: rug_mpi ! rugosity at global grid
204	REAL, DIMENSION(klon_mpi) :: alb_mpi ! albedo at global grid
205
206	CHARACTER(len=20) :: modname='limit_read_mod'
207	CHARACTER(LEN=99) :: abort_message, calendar, str
208
209	! End declaration
210	!****************************************************************************************
211
212	!****************************************************************************************
213	! 0) Initialization
214	!
215	!****************************************************************************************
216	IF (first_call) THEN
217	first_call=.FALSE.
218	! calculate number of time steps for one day
219	lmt_pas = NINT(86400./dtime * 1.0)
220
221	! Allocate module save variables
222	IF ( type_ocean /= 'couple' ) THEN
223	ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
224	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
225	END IF
226
227	IF ( .NOT. ok_veget ) THEN
228	ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
229	IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
230	END IF
231
232	!$OMP MASTER ! Only master thread
233	IF (is_mpi_root) THEN ! Only master processus
234	ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
235	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
236	'Pb d''ouverture du fichier de conditions aux limites',1)
237
238	!--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
239	ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
240	ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
241	IF(ierr==NF90_NOERR.AND.calendar/=calend.AND.prt_level>=1) THEN
242	WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
243	WRITE(lunout,*)' '//TRIM(calend)//' for gcm'
244	WRITE(lunout,*)' '//TRIM(calendar)//' for limit.nc file'
245	END IF
246
247	!--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS
248	IF (grid_type==unstructured) THEN
249	ierr=NF90_INQ_DIMID(nid,"time_year",ndimid)
250	ELSE
251	ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
252	ENDIF
253	ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
254	WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
255	't match year length (',year_len,')'
256	IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
257
258	!--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
259	IF (grid_type==unstructured) THEN
260	ierr=NF90_INQ_DIMID(nid, 'cell', ndimid)
261	ELSE
262	ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
263	ENDIF
264	ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
265	WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
266	') does not match LMDZ klon_glo (',klon_glo,')'
267	IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
268
269	ierr = NF90_CLOSE(nid)
270	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
271	END IF ! is_mpi_root
272	!$OMP END MASTER
273	!$OMP BARRIER
274	END IF
275
276	!****************************************************************************************
277	! 1) Open the file limit.nc if it is the right moment to read, once a day.
278	! The file is read only by the master thread of the master mpi process(is_mpi_root)
279	! Check by the way if the number of records is correct.
280	!
281	!****************************************************************************************
282
283	is_modified = .FALSE.
284	!ym IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN ! time to read
285	! not REALLY PERIODIC
286	IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read
287	! IF (MOD(itime-1, lmt_pas) == 0) THEN ! time to read
288	jour_lu = jour
289	is_modified = .TRUE.
290
291	IF (grid_type==unstructured) THEN
292
293	IF ( type_ocean /= 'couple') THEN
294
295	IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce))
296	IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic))
297	! IF (read_continents .OR. itime == 1) THEN
298	IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter))
299	IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic))
300	! ENDIF
301	ENDIF! type_ocean /= couple
302
303	IF ( type_ocean /= 'couple') THEN
304	IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
305	ENDIF
306
307	IF (.NOT. ok_veget) THEN
308	IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi)
309	IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi)
310	ENDIF
311
312	IF ( type_ocean /= 'couple') THEN
313	CALL Scatter_omp(sst_mpi,sst)
314	CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
315	CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
316	! IF (read_continents .OR. itime == 1) THEN
317	CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter))
318	CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic))
319	! END IF
320	END IF
321
322	IF (.NOT. ok_veget) THEN
323	CALL Scatter_omp(alb_mpi, albedo)
324	CALL Scatter_omp(rug_mpi, rugos)
325	END IF
326
327	ELSE ! grid_type==regular
328
329	!$OMP MASTER ! Only master thread
330	IF (is_mpi_root) THEN ! Only master processus!
331
332	ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
333	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
334	'Pb d''ouverture du fichier de conditions aux limites',1)
335
336	! La tranche de donnees a lire:
337	start(1) = 1
338	start(2) = jour
339	epais(1) = klon_glo
340	epais(2) = 1
341
342
343	!****************************************************************************************
344	! 2) Read fraction if not type_ocean=couple
345	!
346	!****************************************************************************************
347
348	IF ( type_ocean /= 'couple') THEN
349	!
350	! Ocean fraction
351	ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
352	IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
353
354	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
355	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
356	!
357	! Sea-ice fraction
358	ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
359	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
360
361	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
362	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
363
364
365	! Read land and continentals fraction only if asked for
366	IF (read_continents .OR. itime == 1) THEN
367	!
368	! Land fraction
369	ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
370	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
371
372	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
373	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
374	!
375	! Continentale ice fraction
376	ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
377	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
378
379	ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
380	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
381	END IF
382
383	END IF ! type_ocean /= couple
384
385	!****************************************************************************************
386	! 3) Read sea-surface temperature, if not coupled ocean
387	!
388	!****************************************************************************************
389	IF ( type_ocean /= 'couple') THEN
390
391	ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
392	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
393
394	ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
395	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
396
397	END IF
398
399	!****************************************************************************************
400	! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
401	!
402	!****************************************************************************************
403
404	IF (.NOT. ok_veget) THEN
405	!
406	! Read albedo
407	ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
408	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
409
410	ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
411	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
412	!
413	! Read rugosity
414	ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
415	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
416
417	ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
418	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
419
420	END IF
421
422	!****************************************************************************************
423	! 5) Close file and distribuate variables to all processus
424	!
425	!****************************************************************************************
426	ierr = NF90_CLOSE(nid)
427	IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
428	ENDIF ! is_mpi_root
429
430	!$OMP END MASTER
431	!$OMP BARRIER
432
433	IF ( type_ocean /= 'couple') THEN
434	CALL Scatter(sst_glo,sst)
435	CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
436	CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
437	IF (read_continents .OR. itime == 1) THEN
438	CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
439	CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
440	END IF
441	END IF
442
443	IF (.NOT. ok_veget) THEN
444	CALL Scatter(alb_glo, albedo)
445	CALL Scatter(rug_glo, rugos)
446	END IF
447
448	ENDIF ! Grid type
449
450	ENDIF ! time to read
451
452	END SUBROUTINE limit_read_tot
453
454	END MODULE limit_read_mod

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

Download in other formats:

Original Format