Context Navigation

readaerosolstrato2_rrtm.F90 @ 3372

Last change on this file since 3372 was 2744, checked in by oboucher, 8 years ago
Adding threadprivate for mth_pre variable This should fix some reproducibility effects for runs longer than 1 months....
File size: 13.8 KB

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

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

Download in other formats:

Original Format