Context Navigation

readaerosolstrato2_rrtm.F90 @ 3435

Last change on this file since 3435 was 3435, checked in by Laurent Fairhead, 5 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

File size: 15.7 KB

Line
1	!
2	! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
3	!
4
5	SUBROUTINE readaerosolstrato2_rrtm(debut)
6
7	USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
8	nf95_inq_varid, nf95_open
9	USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
10
11	USE phys_cal_mod, ONLY : mth_cur
12	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
13	USE mod_phys_lmdz_mpi_data
14	USE mod_phys_lmdz_omp_data
15	USE mod_phys_lmdz_para
16	USE phys_state_var_mod
17	USE phys_local_var_mod
18	USE aero_mod
19	USE dimphy
20	USE YOERAD, ONLY : NLW
21	USE YOMCST
22	USE xios
23	IMPLICIT NONE
24
25	INCLUDE "clesphys.h"
26
27	CHARACTER (len = 80) :: abort_message
28	CHARACTER (LEN=20) :: modname = 'readaerosolstrato2'
29
30	! Variable input
31	LOGICAL, INTENT(IN) :: debut
32
33	! Variables locales
34	INTEGER n_lat ! number of latitudes in the input data
35	INTEGER n_lon ! number of longitudes
36	INTEGER n_lev ! number of levels in the input data
37	INTEGER n_month ! number of months in the input data
38	INTEGER n_wav ! number of wavelengths in the input data
39	REAL, POINTER:: latitude(:)
40	REAL, POINTER:: time(:)
41	REAL, POINTER:: lev(:)
42	REAL, POINTER:: wav(:)
43	INTEGER i,k,wave,band
44	INTEGER, SAVE :: mth_pre=1
45	!$OMP THREADPRIVATE(mth_pre)
46
47	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: tau_aer_strat
48	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: piz_aer_strat
49	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cg_aer_strat
50	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: taulw_aer_strat
51	!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat)
52
53	! Champs reconstitues
54	REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
55	REAL, ALLOCATABLE:: pizaerstrat(:, :, :, :)
56	REAL, ALLOCATABLE:: cgaerstrat(:, :, :, :)
57	REAL, ALLOCATABLE:: taulwaerstrat(:, :, :, :)
58
59	REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :, :)
60	REAL, ALLOCATABLE:: pizaerstrat_mois(:, :, :, :)
61	REAL, ALLOCATABLE:: cgaerstrat_mois(:, :, :, :)
62	REAL, ALLOCATABLE:: taulwaerstrat_mois(:, :, :, :)
63
64	REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :, :)
65	REAL, ALLOCATABLE:: pizaerstrat_mois_glo(:, :, :)
66	REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
67	REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
68	REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
69	REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
70	REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
71	REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
72
73	! For NetCDF:
74	INTEGER ncid_in ! IDs for input files
75	INTEGER varid, ncerr
76
77	!--------------------------------------------------------
78
79	IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,NSW))
80	IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,NSW))
81	IF (.not.ALLOCATED(cg_aer_strat)) ALLOCATE(cg_aer_strat(klon,klev,NSW))
82
83	IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,NLW))
84
85	!--we only read monthly strat aerosol data
86	IF (debut.OR.mth_cur.NE.mth_pre) THEN
87
88	!--only root reads the data
89	IF (is_mpi_root.AND.is_omp_root) THEN
90
91	!--check mth_cur
92	IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
93	print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
94	ENDIF
95
96	!--initialize n_lon as input data is 2D (lat-alt) only
97	n_lon = nbp_lon
98
99	!--Starts with SW optical properties
100
101	CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
102
103	CALL nf95_inq_varid(ncid_in, "LEV", varid)
104	CALL nf95_gw_var(ncid_in, varid, lev)
105	n_lev = size(lev)
106	IF (n_lev.NE.klev) THEN
107	abort_message='Le nombre de niveaux n est pas egal a klev'
108	CALL abort_physic(modname,abort_message,1)
109	ENDIF
110
111	CALL nf95_inq_varid(ncid_in, "LAT", varid)
112	CALL nf95_gw_var(ncid_in, varid, latitude)
113	n_lat = size(latitude)
114
115	IF (grid_type/=unstructured) THEN
116	IF (n_lat.NE.nbp_lat) THEN
117	print *, 'latitude=', n_lat, nbp_lat
118	abort_message='Le nombre de lat n est pas egal a nbp_lat'
119	CALL abort_physic(modname,abort_message,1)
120	ENDIF
121	ENDIF
122
123	CALL nf95_inq_varid(ncid_in, "TIME", varid)
124	CALL nf95_gw_var(ncid_in, varid, time)
125	n_month = size(time)
126	IF (n_month.NE.12) THEN
127	abort_message='Le nombre de month n est pas egal a 12'
128	CALL abort_physic(modname,abort_message,1)
129	ENDIF
130
131	CALL nf95_inq_varid(ncid_in, "WAV", varid)
132	CALL nf95_gw_var(ncid_in, varid, wav)
133	n_wav = size(wav)
134	print *, 'WAV aerosol strato=', n_wav, wav
135	IF (n_wav.NE.NSW) THEN
136	abort_message='Le nombre de wav n est pas egal a NSW'
137	CALL abort_physic(modname,abort_message,1)
138	ENDIF
139
140	ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
141	ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
142	ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
143
144	!--reading stratospheric aerosol tau per layer
145	CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
146	ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
147	print *,'code erreur readaerosolstrato=', ncerr, varid
148
149	!--reading stratospheric aerosol omega per layer
150	CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
151	ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
152	print *,'code erreur readaerosolstrato=', ncerr, varid
153
154	!--reading stratospheric aerosol g per layer
155	CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
156	ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
157	print *,'code erreur readaerosolstrato sw=', ncerr, varid
158
159	CALL nf95_close(ncid_in)
160
161
162	IF (grid_type/=unstructured) THEN
163	ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
164	ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
165	ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
166
167	ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
168	ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
169	ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
170	!--select the correct month
171	!--and copy into 1st longitude
172	tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
173	pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
174	cgaerstrat_mois(1,:,:,:) = cgaerstrat(:,:,:,mth_cur)
175
176	!--copy longitudes
177	DO i=2, n_lon
178	tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
179	pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
180	cgaerstrat_mois(i,:,:,:) = cgaerstrat_mois(1,:,:,:)
181	ENDDO
182
183	!---reduce to a klon_glo grid
184	DO band=1, NSW
185	CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
186	CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
187	CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
188	ENDDO
189	ENDIF
190	!--Now LW optical properties
191	!
192
193	CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
194
195	CALL nf95_inq_varid(ncid_in, "LEV", varid)
196	CALL nf95_gw_var(ncid_in, varid, lev)
197	n_lev = size(lev)
198	IF (n_lev.NE.klev) THEN
199	abort_message='Le nombre de niveaux n est pas egal a klev'
200	CALL abort_physic(modname,abort_message,1)
201	ENDIF
202
203	CALL nf95_inq_varid(ncid_in, "LAT", varid)
204	CALL nf95_gw_var(ncid_in, varid, latitude)
205	n_lat = size(latitude)
206
207	IF (grid_type/=unstructured) THEN
208	IF (n_lat.NE.nbp_lat) THEN
209	abort_message='Le nombre de lat n est pas egal a nbp_lat'
210	CALL abort_physic(modname,abort_message,1)
211	ENDIF
212	ENDIF
213
214	CALL nf95_inq_varid(ncid_in, "TIME", varid)
215	CALL nf95_gw_var(ncid_in, varid, time)
216	n_month = size(time)
217	IF (n_month.NE.12) THEN
218	abort_message='Le nombre de month n est pas egal a 12'
219	CALL abort_physic(modname,abort_message,1)
220	ENDIF
221
222	CALL nf95_inq_varid(ncid_in, "WAV", varid)
223	CALL nf95_gw_var(ncid_in, varid, wav)
224	n_wav = size(wav)
225	print *, 'WAV aerosol strato=', n_wav, wav
226	IF (n_wav.NE.NLW) THEN
227	abort_message='Le nombre de wav n est pas egal a NLW'
228	CALL abort_physic(modname,abort_message,1)
229	ENDIF
230
231	ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
232
233	!--reading stratospheric aerosol lw tau per layer
234	CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
235	ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
236	print *,'code erreur readaerosolstrato lw=', ncerr, varid
237
238	CALL nf95_close(ncid_in)
239
240	IF (grid_type/=unstructured) THEN
241
242	ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
243	ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
244
245	!--select the correct month
246	!--and copy into 1st longitude
247	taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
248	!--copy longitudes
249	DO i=2, n_lon
250	taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
251	ENDDO
252
253	!---reduce to a klon_glo grid
254	DO band=1, NLW
255	CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
256	ENDDO
257	ENDIF
258
259	ELSE !--proc other than mpi_root and omp_root
260	!--dummy allocation needed for debug mode
261
262	ALLOCATE(tauaerstrat_mois_glo(1,1,1))
263	ALLOCATE(pizaerstrat_mois_glo(1,1,1))
264	ALLOCATE(cgaerstrat_mois_glo(1,1,1))
265	ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
266
267	ALLOCATE(tauaerstrat(0,0,0,12))
268	ALLOCATE(pizaerstrat(0,0,0,12))
269	ALLOCATE(cgaerstrat(0,0,0,12))
270	ALLOCATE(taulwaerstrat(0,0,0,12))
271
272
273	ENDIF !--is_mpi_root and is_omp_root
274
275	!$OMP BARRIER
276
277	!--keep memory of previous month
278	mth_pre=mth_cur
279
280	IF (grid_type==unstructured) THEN
281
282	IF (is_omp_master) THEN
283	ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, NSW))
284	ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, NSW))
285	ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, NSW))
286	ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, NLW))
287
288	CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
289	CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
290	CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
291	CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
292	CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
293	CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
294	CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
295	CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
296	ELSE
297	ALLOCATE(tauaerstrat_mpi(0, 0, 0))
298	ALLOCATE(pizaerstrat_mpi(0, 0, 0))
299	ALLOCATE(cgaerstrat_mpi(0, 0, 0))
300	ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
301	ENDIF
302
303	CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
304	CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
305	CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
306	CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
307	ELSE
308
309	!--scatter on all proc
310	CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
311	CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
312	CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
313	CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
314
315	ENDIF
316
317	IF (is_mpi_root.AND.is_omp_root) THEN
318	DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
319	DEALLOCATE(taulwaerstrat_mois)
320	DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat,taulwaerstrat)
321	ENDIF
322
323
324	!$OMP BARRIER
325
326	ENDIF !--debut ou nouveau mois
327
328	!--total vertical aod at the 5 SW wavelengths
329	!--for now use band 3 AOD into all 5 wavelengths
330	!--it is only a reasonable approximation for 550 nm (wave=2)
331	band=3
332	DO i=1, klon
333	DO k=1, klev
334	IF (stratomask(i,k).GT.0.999999) THEN
335	DO wave=1, nwave_sw
336	tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
337	ENDDO
338	ENDIF
339	ENDDO
340	ENDDO
341
342	!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
343	DO band=1, NSW
344	WHERE (stratomask.GT.0.999999)
345	!--anthropogenic aerosols bands 1 to NSW
346	cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)piz_aero_sw_rrtm(:,:,2,band)tau_aero_sw_rrtm(:,:,2,band) + &
347	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
348	MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
349	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
350	piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
351	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
352	MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
353	tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
354	!--natural aerosols bands 1 to NSW
355	cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)piz_aero_sw_rrtm(:,:,1,band)tau_aero_sw_rrtm(:,:,1,band) + &
356	cg_aer_strat(:,:,band)piz_aer_strat(:,:,band)tau_aer_strat(:,:,band) ) / &
357	MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
358	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
359	piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
360	piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) / &
361	MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
362	tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
363	!--no stratospheric aerosol in index 1 for these tests
364	! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band)
365	! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band)
366	! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band)
367	ENDWHERE
368	ENDDO
369
370	!--total vertical aod at 10 um
371	!--this is approximated from band 7 of RRTM
372	band=7
373	DO i=1, klon
374	DO k=1, klev
375	IF (stratomask(i,k).GT.0.999999) THEN
376	DO wave=1, nwave_lw
377	tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
378	ENDDO
379	ENDIF
380	ENDDO
381	ENDDO
382
383	DO band=1, NLW
384	WHERE (stratomask.GT.0.999999)
385	tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
386	tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
387	!--no stratospheric aerosols in index 1 for these tests
388	! tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band)
389	ENDWHERE
390	ENDDO
391
392	!--default SSA value if there is no aerosol
393	!--to avoid 0 values that seems to cause some problem to RRTM
394	WHERE (tau_aero_sw_rrtm.LT.1.e-14)
395	piz_aero_sw_rrtm = 1.0
396	ENDWHERE
397
398	!--in principle this should not be necessary
399	!--as these variables have min values already but just in case
400	!--put 1e-15 min value to both SW and LW AOD
401	tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
402	tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
403
404	END SUBROUTINE readaerosolstrato2_rrtm

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

Download in other formats:

Original Format