Context Navigation

readaerosol.F90 @ 1179

Last change on this file since 1179 was 1179, checked in by jghattas, 15 years ago

Ajout de l'interpolation vertical pour les nouveaux fichiers de forcage des aerosols. Utilisant les anciennes fichiers de SO4 pas d'interpolation possible. Convergence numerique avec la version precedente en utilisant les anciens fichiers des SO4. aerosol_optic.F90 change du nom pour readaerosol_optic.F90 (lecture d'aerosol + optic) Les fichiers de forcage aerosol doit maintenant avoir le suffix .nc.

Correction des bugs pour inca et certain diagnostiques optionelles de radlwsw.
Ajout de test pour le choix advection schema.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 20.2 KB

Line
1	! $Id: readaerosol.F90 1179 2009-06-11 14:18:47Z jghattas $
2	!
3	MODULE readaerosol_mod
4
5	REAL, SAVE :: not_valid=-333.
6
7	CONTAINS
8
9	SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
10
11	!****************************************************************************************
12	! This routine will read the aersosol from file.
13	!
14	! Read a year data with get_aero_fromfile depending on aer_type :
15	! - actuel : read year 1980
16	! - preind : read natural data
17	! - scenario : read one or two years and do eventually linare time interpolation
18	!
19	! Return pointer, pt_out, to the year read or result from interpolation
20	!****************************************************************************************
21	USE dimphy
22
23	IMPLICIT NONE
24
25	INCLUDE "iniprint.h"
26
27	! Input arguments
28	CHARACTER(len=7), INTENT(IN) :: name_aero
29	CHARACTER(len=*), INTENT(IN) :: type ! correspond to aer_type in clesphys.h
30	INTEGER, INTENT(IN) :: iyr_in
31
32	! Output
33	INTEGER, INTENT(OUT) :: klev_src
34	REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels
35	REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels
36	REAL, POINTER, DIMENSION(:,:,:) :: pt_out ! The massvar distributions, DIMENSION(klon, klev_src, 12)
37	REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf ! Surface pression for 12 months
38	REAL, DIMENSION(klon,12), INTENT(OUT) :: load ! Aerosol mass load in each column for 12 months
39
40	! Local variables
41	CHARACTER(len=4) :: cyear
42	REAL, POINTER, DIMENSION(:,:,:) :: pt_2
43	REAL, DIMENSION(klon,12) :: psurf2, load2
44	REAL :: p0 ! Reference pressure
45	INTEGER :: iyr1, iyr2, klev_src2
46	INTEGER :: it, k, i
47	LOGICAL, PARAMETER :: lonlyone=.FALSE.
48
49	!****************************************************************************************
50	! Read data depending on aer_type
51	!
52	!****************************************************************************************
53
54	IF (type == 'actuel') THEN
55	! Read and return data for year 1980
56	!****************************************************************************************
57	cyear='1980'
58	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
59	! pt_out has dimensions (klon, klev_src, 12)
60	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
61
62
63	ELSE IF (type == 'preind') THEN
64	! Read and return data from file with suffix .nat
65	!****************************************************************************************
66	cyear='.nat'
67	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
68	! pt_out has dimensions (klon, klev_src, 12)
69	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
70
71	ELSE IF (type == 'scenario') THEN
72	! Read data depending on actual year and interpolate if necessary
73	!****************************************************************************************
74	IF (iyr_in .LT. 1850) THEN
75	cyear='.nat'
76	WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,' ',cyear
77	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
78	! pt_out has dimensions (klon, klev_src, 12)
79	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
80
81	ELSE IF (iyr_in .GE. 2100) THEN
82	cyear='2100'
83	WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear
84	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
85	! pt_out has dimensions (klon, klev_src, 12)
86	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
87
88	ELSE
89	! Read data from 2 decades and interpolate to actual year
90	! a) from actual 10-yr-period
91	IF (iyr_in.LT.1900) THEN
92	iyr1 = 1850
93	iyr2 = 1900
94	ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN
95	iyr1 = 1900
96	iyr2 = 1920
97	ELSE
98	iyr1 = INT(iyr_in/10)*10
99	iyr2 = INT(1+iyr_in/10)*10
100	ENDIF
101
102	WRITE(cyear,'(I4)') iyr1
103	WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear
104	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
105	! pt_out has dimensions (klon, klev_src, 12)
106	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
107
108	! If to read two decades:
109	IF (.NOT.lonlyone) THEN
110
111	! b) from the next following one
112	WRITE(cyear,'(I4)') iyr2
113	WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,' ',cyear
114
115	NULLIFY(pt_2)
116	! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
117	! pt_2 has dimensions (klon, klev_src, 12)
118	CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
119	! Test for same number of vertical levels
120	IF (klev_src /= klev_src2) THEN
121	WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded'
122	CALL abort_gcm('readaersosol','Error in number of vertical levels',1)
123	END IF
124
125	! Linare interpolate to the actual year:
126	DO it=1,12
127	DO k=1,klev_src
128	DO i = 1, klon
129	pt_out(i,k,it) = &
130	pt_out(i,k,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
131	(pt_out(i,k,it) - pt_2(i,k,it))
132	END DO
133	END DO
134
135	DO i = 1, klon
136	psurf(i,it) = &
137	psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
138	(psurf(i,it) - psurf2(i,it))
139
140	load(i,it) = &
141	load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
142	(load(i,it) - load2(i,it))
143	END DO
144	END DO
145
146	! Deallocate pt_2 no more needed
147	DEALLOCATE(pt_2)
148
149	END IF ! lonlyone
150	END IF ! iyr_in .LT. 1850
151
152	ELSE
153	WRITE(lunout,*)'This option is not implemented : aer_type = ', type
154	CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
155	END IF ! type
156
157
158	END SUBROUTINE readaerosol
159
160
161	SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
162	!****************************************************************************************
163	! Read 12 month aerosol from file and distribute to local process on physical grid.
164	! Vertical levels, klev_src, may differ from model levels if new file format.
165	!
166	! For mpi_root and master thread :
167	! 1) Open file
168	! 2) Find vertical dimension klev_src
169	! 3) Read field month by month
170	! 4) Close file
171	! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo)
172	! - Also the levels and the latitudes have to be inversed
173	!
174	! For all processes and threads :
175	! 6) Scatter global field(klon_glo) to local process domain(klon)
176	! 7) Test for negative values
177	!****************************************************************************************
178
179	USE netcdf
180	USE dimphy
181	USE mod_grid_phy_lmdz
182	USE mod_phys_lmdz_para
183
184	IMPLICIT NONE
185
186	INCLUDE "dimensions.h"
187	INCLUDE "iniprint.h"
188
189	! Input argumets
190	CHARACTER(len=7), INTENT(IN) :: varname
191	CHARACTER(len=4), INTENT(IN) :: cyr
192
193	! Output arguments
194	INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file
195	REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels
196	REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels
197	REAL :: p0 ! Reference pressure value
198	REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src
199	REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months
200	REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column
201
202	! Local variables
203	CHARACTER(len=30) :: fname
204	CHARACTER(len=30) :: cvar
205	INTEGER :: ncid, dimid, varid
206	INTEGER :: imth, i, j, k, ierr
207	REAL :: npole, spole
208	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth
209	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month
210	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12)
211	REAL, ALLOCATABLE, DIMENSION(:) :: varktmp
212
213	REAL, DIMENSION(iim,jjm+1,12) :: psurf_glo2D ! Surface pression for 12 months on dynamics global grid
214	REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid
215	REAL, DIMENSION(iim,jjm+1,12) :: load_glo2D ! Load for 12 months on dynamics global grid
216	REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid
217	REAL, DIMENSION(iim,jjm+1) :: vartmp
218	LOGICAL :: new_file
219
220
221	! Deallocate pointers
222	IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
223	IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b)
224
225	!$OMP MASTER
226	IF (is_mpi_root) THEN
227
228	! 1) Open file
229	!****************************************************************************************
230	fname = TRIM(varname)//'.run'//cyr//'.nc'
231
232	WRITE(lunout,*) 'reading ', TRIM(fname)
233	CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
234
235	! 2) Check if old or new file is avalabale.
236	! New type of file should contain the dimension 'lev'
237	! Old type of file should contain the dimension 'PRESNIVS'
238	!****************************************************************************************
239	ierr = nf90_inq_dimid(ncid, 'lev', dimid)
240	IF (ierr /= NF90_NOERR) THEN
241	! Coordinate axe lev not found. Check for presnivs.
242	ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
243	IF (ierr /= NF90_NOERR) THEN
244	! Dimension PRESNIVS not found either
245	CALL abort_gcm('get_aero_fromfile', 'dimension lev or presnivs not in file',1)
246	ELSE
247	! Old file found
248	new_file=.FALSE.
249	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
250	END IF
251	ELSE
252	! New file found
253	new_file=.TRUE.
254	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
255	END IF
256
257	! 2) Find vertical dimension klev_src
258	!****************************************************************************************
259	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) )
260
261	! Allocate variables depending on the number of vertical levels
262	ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
263	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
264
265	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
266	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 2',1)
267
268	! 3) Read all variables from file
269	! There is 2 options for the file structure :
270	! new_file=TRUE : read varyear, ps, pt_ap and pt_b
271	! new_file=FALSE : read varyear month by month
272	!****************************************************************************************
273
274	IF (new_file) THEN
275
276	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
277	!****************************************************************************************
278	! Get variable id
279	CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid) )
280
281	! Get the variable
282	CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
283
284	! ++) Read surface pression, 12 month in one variable
285	!****************************************************************************************
286	! Get variable id
287	CALL check_err( nf90_inq_varid(ncid, "ps", varid) )
288	! Get the variable
289	CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) )
290
291	! ++) Read mass load, 12 month in one variable
292	!****************************************************************************************
293	! Get variable id
294	CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) )
295	! Get the variable
296	CALL check_err( nf90_get_var(ncid, varid, load_glo2D) )
297
298	! ++) Read ap
299	!****************************************************************************************
300	! Get variable id
301	CALL check_err( nf90_inq_varid(ncid, "ap", varid) )
302	! Get the variable
303	CALL check_err( nf90_get_var(ncid, varid, pt_ap) )
304
305	! ++) Read b
306	!****************************************************************************************
307	! Get variable id
308	CALL check_err( nf90_inq_varid(ncid, "b", varid) )
309	! Get the variable
310	CALL check_err( nf90_get_var(ncid, varid, pt_b) )
311
312	! ++) Read p0 : reference pressure
313	!****************************************************************************************
314	! Get variable id
315	CALL check_err( nf90_inq_varid(ncid, "p0", varid) )
316	! Get the variable
317	CALL check_err( nf90_get_var(ncid, varid, p0) )
318
319
320	ELSE ! old file
321
322	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
323	!****************************************************************************************
324	DO imth=1, 12
325	IF (imth.EQ.1) THEN
326	cvar=TRIM(varname)//'JAN'
327	ELSE IF (imth.EQ.2) THEN
328	cvar=TRIM(varname)//'FEB'
329	ELSE IF (imth.EQ.3) THEN
330	cvar=TRIM(varname)//'MAR'
331	ELSE IF (imth.EQ.4) THEN
332	cvar=TRIM(varname)//'APR'
333	ELSE IF (imth.EQ.5) THEN
334	cvar=TRIM(varname)//'MAY'
335	ELSE IF (imth.EQ.6) THEN
336	cvar=TRIM(varname)//'JUN'
337	ELSE IF (imth.EQ.7) THEN
338	cvar=TRIM(varname)//'JUL'
339	ELSE IF (imth.EQ.8) THEN
340	cvar=TRIM(varname)//'AUG'
341	ELSE IF (imth.EQ.9) THEN
342	cvar=TRIM(varname)//'SEP'
343	ELSE IF (imth.EQ.10) THEN
344	cvar=TRIM(varname)//'OCT'
345	ELSE IF (imth.EQ.11) THEN
346	cvar=TRIM(varname)//'NOV'
347	ELSE IF (imth.EQ.12) THEN
348	cvar=TRIM(varname)//'DEC'
349	END IF
350
351	! Get variable id
352	CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid) )
353
354	! Get the variable
355	CALL check_err( nf90_get_var(ncid, varid, varmth) )
356
357	! Store in variable for the whole year
358	varyear(:,:,:,imth)=varmth(:,:,:)
359
360	END DO
361
362	! Putting dummy
363	psurf_glo2D(:,:,:) = not_valid
364	load_glo2D(:,:,:) = not_valid
365	pt_ap(:) = not_valid
366	pt_b(:) = not_valid
367
368	END IF
369
370	! 4) Close file
371	!****************************************************************************************
372	CALL check_err( nf90_close(ncid) )
373
374
375	! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo)
376	!****************************************************************************************
377	! Test if vertical levels have to be inversed
378
379	IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
380	WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inversed'
381	WRITE(lunout,*) 'before pt_ap = ', pt_ap
382	WRITE(lunout,*) 'before pt_b = ', pt_b
383
384	! Inverse vertical levels for varyear
385	DO imth=1, 12
386	varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
387	DO k=1, klev_src
388	DO j=1, jjm+1
389	DO i=1,iim
390	varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
391	END DO
392	END DO
393	END DO
394	END DO
395
396	! Inverte vertical axes for pt_ap and pt_b
397	varktmp(:) = pt_ap(:)
398	DO k=1, klev_src
399	pt_ap(k) = varktmp(klev_src+1-k)
400	END DO
401
402	varktmp(:) = pt_b(:)
403	DO k=1, klev_src
404	pt_b(k) = varktmp(klev_src+1-k)
405	END DO
406	WRITE(lunout,*) 'after pt_ap = ', pt_ap
407	WRITE(lunout,*) 'after pt_b = ', pt_b
408
409	ELSE
410	WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
411	WRITE(lunout,*) 'pt_ap = ', pt_ap
412	WRITE(lunout,*) 'pt_b = ', pt_b
413	END IF
414
415	! - Also the latitudes have to be inversed
416	DO imth=1, 12
417	! Inverse latitudes
418	varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
419	DO k=1,klev_src
420	DO j=1,jjm+1
421	DO i=1,iim
422	varyear(i,j,k,imth) = varmth(i,jjm+1+1-j,k)
423	END DO
424	END DO
425	END DO
426
427	! Inverse latitudes for surface pressure
428	vartmp(:,:) = psurf_glo2D(:,:,imth)
429	DO j=1, jjm+1
430	DO i=1,iim
431	psurf_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
432	END DO
433	END DO
434
435	! Inverse latitudes for the load
436	vartmp(:,:) = load_glo2D(:,:,imth)
437	DO j=1, jjm+1
438	DO i=1,iim
439	load_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
440	END DO
441	END DO
442
443	! Do zonal mead at poles and distribut at whole first and last latitude
444	DO k=1, klev_src
445	npole=0. ! North pole, j=1
446	spole=0. ! South pole, j=jjm+1
447	DO i=1,iim
448	npole = npole + varyear(i,1,k,imth)
449	spole = spole + varyear(i,jjm+1,k,imth)
450	END DO
451	npole = npole/FLOAT(iim)
452	spole = spole/FLOAT(iim)
453	varyear(:,1, k,imth) = npole
454	varyear(:,jjm+1,k,imth) = spole
455	END DO
456	END DO ! imth
457
458	ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
459	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
460
461	! Transform from 2D to 1D field
462	CALL grid2Dto1D_glo(varyear,varyear_glo1D)
463	CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
464	CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
465
466	END IF ! is_mpi_root
467	!$OMP END MASTER BARRIER
468
469	! 6) Distribute to all processes
470	! Scatter global field(klon_glo) to local process domain(klon)
471	! and distribute klev_src to all processes
472	!****************************************************************************************
473
474	! Distribute klev_src
475	CALL bcast(klev_src)
476
477	! Allocate and distribute pt_ap and pt_b
478	IF (.NOT. ASSOCIATED(pt_ap)) THEN ! if pt_ap is allocated also pt_b is allocated
479	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
480	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 4',1)
481	END IF
482	CALL bcast(pt_ap)
483	CALL bcast(pt_b)
484
485	! Allocate space for output pointer variable at local process
486	IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
487	ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
488	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1)
489
490	! Scatter global field to local domain at local process
491	CALL scatter(varyear_glo1D, pt_year)
492	CALL scatter(psurf_glo1D, psurf_out)
493	CALL scatter(load_glo1D, load_out)
494
495	! 7) Test for negative values
496	!****************************************************************************************
497	IF (MINVAL(pt_year) < 0.) THEN
498	WRITE(lunout,*) 'Warning! Negative values read from file :', fname
499	END IF
500
501	END SUBROUTINE get_aero_fromfile
502
503
504	SUBROUTINE check_err(status)
505	USE netcdf
506	IMPLICIT NONE
507
508	INCLUDE "iniprint.h"
509	INTEGER, INTENT (IN) :: status
510
511	IF (status /= NF90_NOERR) THEN
512	WRITE(lunout,*) 'Error in get_aero_fromfile ',status
513	CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
514	END IF
515
516	END SUBROUTINE check_err
517
518
519	END MODULE readaerosol_mod

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

Download in other formats:

Original Format