Context Navigation

readaerosolstrato2_rrtm.F90 @ 5429

Last change on this file since 5429 was 3427, checked in by oboucher, 6 years ago
Correcting a bug for strat aerosol forcing diagnostics
File size: 15.5 KB

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

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

Download in other formats:

Original Format