Context Navigation

readaerosolstrato1_rrtm.F90 @ 3435

Last change on this file since 3435 was 3435, checked in by Laurent Fairhead, 6 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: 10.2 KB

Line
1	!
2	! $Id: readaerosolstrato1_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
3	!
4
5	SUBROUTINE readaerosolstrato1_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_para
14	USE phys_state_var_mod
15	USE phys_local_var_mod
16	USE aero_mod
17	USE dimphy
18	USE YOERAD, ONLY : NLW
19	USE YOMCST
20	USE xios
21
22	IMPLICIT NONE
23
24	! Variable input
25	LOGICAL debut
26
27	! Variables locales
28	INTEGER n_lat ! number of latitudes in the input data
29	INTEGER n_lon ! number of longitudes in the input data
30	INTEGER n_lev ! number of levels in the input data
31	INTEGER n_month ! number of months in the input data
32	REAL, POINTER:: latitude(:)
33	REAL, POINTER:: longitude(:)
34	REAL, POINTER:: time(:)
35	REAL, POINTER:: lev(:)
36	INTEGER k, band, wave, i
37	INTEGER, SAVE :: mth_pre=1
38	!$OMP THREADPRIVATE(mth_pre)
39
40	REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tau_aer_strat
41	!$OMP THREADPRIVATE(tau_aer_strat)
42
43	! Champs reconstitues
44	REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
45	REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :)
46	REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :)
47	REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :)
48
49	! For NetCDF:
50	INTEGER ncid_in ! IDs for input files
51	INTEGER varid, ncerr
52
53	! Stratospheric aerosols optical properties
54	! alpha_sw_strat over the 6 bands is normalised by the 550 nm extinction coefficient
55	REAL, DIMENSION(nbands_sw_rrtm) :: alpha_sw_strat, piz_sw_strat, cg_sw_strat
56	DATA alpha_sw_strat/0.8545564, 0.8451642, 0.9821724, 0.8145110, 0.3073565, 7.7966176E-02/
57	DATA cg_sw_strat /0.6997170, 0.6810035, 0.7403592, 0.7562674, 0.6676504, 0.3478689/
58	DATA piz_sw_strat /0.9999998, 0.9999998, 1.000000000, 0.9999958, 0.9977155, 0.4510679/
59	!
60	!--diagnostics AOD in the SW
61	! alpha_sw_strat_wave is not normalised by the 550 nm extinction coefficient
62	REAL, DIMENSION(nwave_sw) :: alpha_sw_strat_wave
63	DATA alpha_sw_strat_wave/3.708007,4.125824,4.136584,3.887478,3.507738/
64	!
65	!--diagnostics AOD in the LW at 10 um (not normalised by the 550 nm ext coefficient
66	REAL :: alpha_lw_strat_wave(nwave_lw)
67	DATA alpha_lw_strat_wave/0.2746812/
68	!
69	REAL, DIMENSION(nbands_lw_rrtm) :: alpha_lw_abs_rrtm
70	DATA alpha_lw_abs_rrtm/ 8.8340312E-02, 6.9856711E-02, 6.2652975E-02, 5.7188231E-02, &
71	6.3157059E-02, 5.5072524E-02, 5.0571125E-02, 0.1349073, &
72	0.1381676, 9.6506312E-02, 5.1312990E-02, 2.4256418E-02, &
73	2.7191756E-02, 3.3862915E-02, 1.6132960E-02, 1.4275438E-02/ ! calculated with Mie_SW_LW_RRTM_V2.4 (bimodal, corrected)
74	! for r_0=/0.13E-6, 0.41E-6/ m, sigma_g=/1.26, 1.30/
75	! order: increasing wavelength!
76	!--------------------------------------------------------
77
78	IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev))
79
80	!--we only read monthly strat aerosol data
81	IF (debut.OR.mth_cur.NE.mth_pre) THEN
82
83	!--only root reads the data
84	IF (is_mpi_root.AND.is_omp_root) THEN
85
86	IF (nbands_sw_rrtm.NE.6) THEN
87	print *,'nbands_sw_rrtm doit etre egal a 6 dans readaerosolstrat_rrtm'
88	STOP
89	ENDIF
90
91	CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
92
93	CALL nf95_inq_varid(ncid_in, "LEV", varid)
94	CALL nf95_gw_var(ncid_in, varid, lev)
95	n_lev = size(lev)
96	IF (n_lev.NE.klev) THEN
97	print *,'Le nombre de niveaux n est pas egal a klev'
98	STOP
99	ENDIF
100
101	CALL nf95_inq_varid(ncid_in, "LAT", varid)
102	CALL nf95_gw_var(ncid_in, varid, latitude)
103	n_lat = size(latitude)
104	print *, 'LAT aerosol strato=', n_lat, latitude
105
106	IF (grid_type/=unstructured) THEN
107	IF (n_lat.NE.nbp_lat) THEN
108	print *,'Le nombre de lat n est pas egal a nbp_lat'
109	STOP
110	ENDIF
111	ENDIF
112
113	CALL nf95_inq_varid(ncid_in, "LON", varid)
114	CALL nf95_gw_var(ncid_in, varid, longitude)
115	n_lon = size(longitude)
116	print *, 'LON aerosol strato=', n_lon, longitude
117
118	IF (grid_type/=unstructured) THEN
119	IF (n_lon.NE.nbp_lon) THEN
120	print *,'Le nombre de lon n est pas egal a nbp_lon'
121	STOP
122	ENDIF
123	ENDIF
124
125
126	CALL nf95_inq_varid(ncid_in, "TIME", varid)
127	CALL nf95_gw_var(ncid_in, varid, time)
128	n_month = size(time)
129	print *, 'TIME aerosol strato=', n_month, time
130	IF (n_month.NE.12) THEN
131	print *,'Le nombre de month n est pas egal a 12'
132	STOP
133	ENDIF
134
135	ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
136	ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
137	ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
138
139	!--reading stratospheric AOD at 550 nm
140	CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
141	ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
142	print *,'code erreur readaerosolstrato=', ncerr, varid
143
144	CALL nf95_close(ncid_in)
145
146	!---select the correct month
147	IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
148	print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
149	ENDIF
150	tauaerstrat_mois(:,:,:) = tauaerstrat(:,:,:,mth_cur)
151
152	!---reduce to a klon_glo grid
153	CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
154
155	ELSE
156	ALLOCATE(tauaerstrat_mois(0,0,0))
157	ENDIF !--is_mpi_root and is_omp_root
158
159	!$OMP BARRIER
160
161	!--keep memory of previous month
162	mth_pre=mth_cur
163
164	!--scatter on all proc
165
166	IF (grid_type==unstructured) THEN
167	IF (is_omp_master) THEN
168	ALLOCATE(tauaerstrat_mpi(klon_mpi,klev))
169	CALL xios_send_field("taustrat_in",tauaerstrat_mois)
170	CALL xios_recv_field("taustrat_out",tauaerstrat_mpi)
171	ELSE
172	ALLOCATE(tauaerstrat_mpi(0,0))
173	ENDIF
174	CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
175	ELSE
176	CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
177	ENDIF
178
179	IF (is_mpi_root.AND.is_omp_root) THEN
180	!
181	DEALLOCATE(tauaerstrat)
182	DEALLOCATE(tauaerstrat_mois)
183	DEALLOCATE(tauaerstrat_mois_glo)
184	!
185	ENDIF !--is_mpi_root and is_omp_root
186
187	!$OMP BARRIER
188
189	ENDIF !--debut ou nouveau mois
190
191	!--total vertical aod at the 5 SW wavelengths
192	DO wave=1, nwave_sw
193	DO k=1, klev
194	tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+ &
195	tau_aer_strat(:,k)*alpha_sw_strat_wave(wave)/alpha_sw_strat_wave(2)
196	ENDDO
197	ENDDO
198
199	!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
200	DO band=1, nbands_sw_rrtm
201	!--anthropogenic aerosols bands 1 to nbands_sw_rrtm
202	cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)piz_aero_sw_rrtm(:,:,2,band)tau_aero_sw_rrtm(:,:,2,band) + &
203	cg_sw_strat(band)piz_sw_strat(band)alpha_sw_strat(band)*tau_aer_strat(:,:) ) / &
204	MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
205	piz_sw_strat(band)alpha_sw_strat(band)tau_aer_strat(:,:), 1.e-15 )
206	piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
207	piz_sw_strat(band)alpha_sw_strat(band)tau_aer_strat(:,:) ) / &
208	MAX( tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
209	tau_aero_sw_rrtm(:,:,2,band) = tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
210	!--natural aerosols bands 1 to nbands_sw_rrtm
211	cg_aero_sw_rrtm(:,:,1,band) = ( cg_aero_sw_rrtm(:,:,1,band)piz_aero_sw_rrtm(:,:,1,band)tau_aero_sw_rrtm(:,:,1,band) + &
212	cg_sw_strat(band)piz_sw_strat(band)alpha_sw_strat(band)*tau_aer_strat(:,:) ) / &
213	MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
214	piz_sw_strat(band)alpha_sw_strat(band)tau_aer_strat(:,:), 1.e-15 )
215	piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
216	piz_sw_strat(band)alpha_sw_strat(band)tau_aer_strat(:,:) ) / &
217	MAX( tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:),1.e-15 )
218	tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
219	!--no stratospheric aerosol in index 1 for these tests
220	! cg_aero_sw_rrtm(:,:,1,band) = cg_aero_sw_rrtm(:,:,1,band)
221	! piz_aero_sw_rrtm(:,:,1,band) = piz_aero_sw_rrtm(:,:,1,band)
222	! tau_aero_sw_rrtm(:,:,1,band) = tau_aero_sw_rrtm(:,:,1,band)
223	ENDDO
224
225	!--stratospheric AOD in LW
226	IF (nbands_lw_rrtm .NE. NLW) then
227	print*, 'different values for NLW (=',NLW,') and nbands_lw_rrtm (=', nbands_lw_rrtm, ')'
228	STOP
229	ENDIF
230
231	!--total vertical aod at the 1 LW wavelength
232	DO wave=1, nwave_lw
233	DO k=1, klev
234	tausum_aero(:,nwave_sw+wave,id_STRAT_phy)=tausum_aero(:,nwave_sw+wave,id_STRAT_phy)+ &
235	tau_aer_strat(:,k)*alpha_lw_strat_wave(wave)/alpha_sw_strat_wave(2)
236	ENDDO
237	ENDDO
238
239	DO band=1, nbands_lw_rrtm
240	tau_aero_lw_rrtm(:,:,2,band) = tau_aero_lw_rrtm(:,:,2,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)
241	tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)
242	!--no stratospheric aerosols in index 1 for these tests
243	! tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,1,band)
244	ENDDO
245
246	!--default SSA value if there is no aerosol
247	!--to avoid 0 values that seems to cause some problem to RRTM
248	WHERE (tau_aero_sw_rrtm.LT.1.e-14)
249	piz_aero_sw_rrtm = 1.0
250	ENDWHERE
251
252	!--in principle this should not be necessary
253	!--as these variables have min values already but just in case
254	!--put 1e-15 min value to both SW and LW AOD
255	tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
256	tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
257
258	END SUBROUTINE readaerosolstrato1_rrtm

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

Download in other formats:

Original Format