Context Navigation

cosp.F90 @ 5308

Last change on this file since 5308 was 1414, checked in by idelkadi, 14 years ago
Passage a la version cosp.v1.3 pour le Lidar et ISCCP Corrections de bugs pour ISCCP et optimisation pour le Lidar
File size: 26.4 KB

Rev	Line
[1262]	1	! (c) British Crown Copyright 2008, the Met Office.
	2	! All rights reserved.
	3	!
	4	! Redistribution and use in source and binary forms, with or without modification, are permitted
	5	! provided that the following conditions are met:
	6	!
	7	! * Redistributions of source code must retain the above copyright notice, this list
	8	! of conditions and the following disclaimer.
	9	! * Redistributions in binary form must reproduce the above copyright notice, this list
	10	! of conditions and the following disclaimer in the documentation and/or other materials
	11	! provided with the distribution.
	12	! * Neither the name of the Met Office nor the names of its contributors may be used
	13	! to endorse or promote products derived from this software without specific prior written
	14	! permission.
	15	!
	16	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
	17	! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
	18	! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
	19	! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	20	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	21	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
	22	! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
	23	! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	24
[1414]	25	!!#include "cosp_defs.h"
[1262]	26	MODULE MOD_COSP
	27	USE MOD_COSP_TYPES
	28	USE MOD_COSP_SIMULATOR
[1327]	29	USE mod_phys_lmdz_para
	30	USE mod_grid_phy_lmdz
[1262]	31	IMPLICIT NONE
	32
	33	CONTAINS
	34
	35
	36	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	37	!--------------------- SUBROUTINE COSP ---------------------------
	38	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	39	SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
	40
	41	! Arguments
	42	integer,intent(in) :: overlap ! overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
	43	integer,intent(in) :: Ncolumns
	44	type(cosp_config),intent(in) :: cfg ! Configuration options
	45	type(cosp_vgrid),intent(in) :: vgrid ! Information on vertical grid of stats
	46	type(cosp_gridbox),intent(inout) :: gbx
	47	type(cosp_subgrid),intent(inout) :: sgx ! Subgrid info
	48	type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
	49	type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
	50	type(cosp_isccp),intent(inout) :: isccp ! Output from ISCCP simulator
	51	type(cosp_misr),intent(inout) :: misr ! Output from MISR simulator
	52	type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
	53	type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
	54
	55	! Local variables
	56	integer :: Npoints ! Number of gridpoints
	57	integer :: Nlevels ! Number of levels
	58	integer :: Nhydro ! Number of hydrometeors
	59	integer :: Niter ! Number of calls to cosp_simulator
	60	integer :: i_first,i_last ! First and last gridbox to be processed in each iteration
	61	integer :: i,j,k,Ni
	62	integer,dimension(2) :: ix,iy
	63	logical :: reff_zero
	64	real :: minv,maxv
	65	real :: maxp,minp
[1327]	66	integer,dimension(:),save, allocatable :: & ! Dimensions nPoints
[1262]	67	seed ! It is recommended that the seed is set to a different value for each model
	68	! gridbox it is called on, as it is possible that the choice of the same
	69	! seed value every time may introduce some statistical bias in the results,
	70	! particularly for low values of NCOL.
[1327]	71	!$OMP THREADPRIVATE(seed)
	72	real,dimension(:),allocatable :: rseed ! It is recommended that the seed is set to a different value for each model
[1262]	73	! Types used in one iteration
	74	type(cosp_gridbox) :: gbx_it
	75	type(cosp_subgrid) :: sgx_it
	76	type(cosp_vgrid) :: vgrid_it
	77	type(cosp_sgradar) :: sgradar_it
	78	type(cosp_sglidar) :: sglidar_it
	79	type(cosp_isccp) :: isccp_it
	80	type(cosp_misr) :: misr_it
	81	type(cosp_radarstats) :: stradar_it
	82	type(cosp_lidarstats) :: stlidar_it
	83
[1327]	84	logical,save :: first_cosp=.TRUE.
	85	!$OMP THREADPRIVATE(first_cosp)
	86
[1262]	87	!++++++++++ Dimensions ++++++++++++
	88	Npoints = gbx%Npoints
	89	Nlevels = gbx%Nlevels
	90	Nhydro = gbx%Nhydro
	91
	92	!++++++++++ Apply sanity checks to inputs ++++++++++
	93	! call cosp_check_input('longitude',gbx%longitude,min_val=0.0,max_val=360.0)
	94	call cosp_check_input('longitude',gbx%longitude,min_val=-180.0,max_val=180.0)
	95	call cosp_check_input('latitude',gbx%latitude,min_val=-90.0,max_val=90.0)
	96	call cosp_check_input('dlev',gbx%dlev,min_val=0.0)
	97	call cosp_check_input('p',gbx%p,min_val=0.0)
	98	call cosp_check_input('ph',gbx%ph,min_val=0.0)
	99	call cosp_check_input('T',gbx%T,min_val=0.0)
	100	call cosp_check_input('q',gbx%q,min_val=0.0)
	101	call cosp_check_input('sh',gbx%sh,min_val=0.0)
	102	call cosp_check_input('dtau_s',gbx%dtau_s,min_val=0.0)
	103	call cosp_check_input('dtau_c',gbx%dtau_c,min_val=0.0)
	104	call cosp_check_input('dem_s',gbx%dem_s,min_val=0.0,max_val=1.0)
	105	call cosp_check_input('dem_c',gbx%dem_c,min_val=0.0,max_val=1.0)
	106	! Point information (Npoints)
	107	call cosp_check_input('land',gbx%land,min_val=0.0,max_val=1.0)
	108	call cosp_check_input('psfc',gbx%psfc,min_val=0.0)
	109	call cosp_check_input('sunlit',gbx%sunlit,min_val=0.0,max_val=1.0)
	110	call cosp_check_input('skt',gbx%skt,min_val=0.0)
	111	! TOTAL and CONV cloud fraction for SCOPS
	112	call cosp_check_input('tca',gbx%tca,min_val=0.0,max_val=1.0)
	113	call cosp_check_input('cca',gbx%cca,min_val=0.0,max_val=1.0)
	114	! Precipitation fluxes on model levels
	115	call cosp_check_input('rain_ls',gbx%rain_ls,min_val=0.0)
	116	call cosp_check_input('rain_cv',gbx%rain_cv,min_val=0.0)
	117	call cosp_check_input('snow_ls',gbx%snow_ls,min_val=0.0)
	118	call cosp_check_input('snow_cv',gbx%snow_cv,min_val=0.0)
	119	call cosp_check_input('grpl_ls',gbx%grpl_ls,min_val=0.0)
	120	! Hydrometeors concentration and distribution parameters
	121	call cosp_check_input('mr_hydro',gbx%mr_hydro,min_val=0.0)
	122	! Effective radius [m]. (Npoints,Nlevels,Nhydro)
	123	call cosp_check_input('Reff',gbx%Reff,min_val=0.0)
	124	reff_zero=.true.
	125	if (any(gbx%Reff > 1.e-8)) then
	126	reff_zero=.false.
	127	! reff_zero == .false.
	128	! and gbx%use_reff == .true. Reff use in radar and lidar
	129	! and reff_zero == .false. Reff use in lidar and set to 0 for radar
	130	endif
[1414]	131	! if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
	132	! print *, '---------- COSP ERROR ------------'
	133	! print *, ''
	134	! print *, 'use_reff==.true. but Reff is always zero'
	135	! print *, ''
	136	! print *, '----------------------------------'
	137	! stop
	138	! endif
[1262]	139	if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
	140	gbx%Reff = DEFAULT_LIDAR_REFF
	141	print *, '---------- COSP WARNING ------------'
	142	print *, ''
	143	print *, 'Using default Reff in lidar simulations'
	144	print *, ''
	145	print *, '----------------------------------'
	146	endif
	147
	148	! Aerosols concentration and distribution parameters
	149	call cosp_check_input('conc_aero',gbx%conc_aero,min_val=0.0)
	150	! Checks for CRM mode
	151	if (Ncolumns == 1) then
	152	if (gbx%use_precipitation_fluxes) then
	153	print *, '---------- COSP ERROR ------------'
	154	print *, ''
	155	print *, 'Use of precipitation fluxes not supported in CRM mode (Ncolumns=1)'
	156	print *, ''
	157	print *, '----------------------------------'
	158	stop
	159	endif
	160	if ((maxval(gbx%dtau_c) > 0.0).or.(maxval(gbx%dem_c) > 0.0)) then
	161	print *, '---------- COSP ERROR ------------'
	162	print *, ''
	163	print *, ' dtau_c > 0.0 or dem_c > 0.0. In CRM mode (Ncolumns=1), '
	164	print *, ' the optical depth (emmisivity) of all clouds must be '
	165	print *, ' passed through dtau_s (dem_s)'
	166	print *, ''
	167	print *, '----------------------------------'
	168	stop
	169	endif
	170	endif
	171
[1327]	172	if (first_cosp) then
[1262]	173	! We base the seed in the decimal part of the surface pressure.
[1327]	174	allocate(seed(Npoints))
	175
	176	allocate(rseed(klon_glo))
	177	CALL gather(gbx%psfc,rseed)
	178	call bcast(rseed)
	179	! seed = int(gbx%psfc) ! This is to avoid division by zero when Npoints = 1
[1262]	180	! Roj Oct/2008 ... Note: seed value of 0 caused me some problems + I want to
	181	! randomize for each call to COSP even when Npoints ==1
[1327]	182	minp = minval(rseed)
	183	maxp = maxval(rseed)
	184
	185	if (Npoints .gt. 1) THEN
	186	seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
	187	else
	188	seed=int(gbx%psfc-minp)
	189	endif
[1262]	190
[1327]	191	deallocate(rseed)
	192	first_cosp=.false.
	193	endif
	194
[1262]	195	if (gbx%Npoints_it >= gbx%Npoints) then ! One iteration gbx%Npoints
	196	call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
	197	else ! Several iterations to save memory
	198	Niter = gbx%Npoints/gbx%Npoints_it ! Integer division
	199	if (Niter*gbx%Npoints_it < gbx%Npoints) Niter = Niter + 1
	200	do i=1,Niter
	201	i_first = (i-1)*gbx%Npoints_it + 1
	202	i_last = i_first + gbx%Npoints_it - 1
	203	i_last = min(i_last,gbx%Npoints)
	204	Ni = i_last - i_first + 1
	205	if (i == 1) then
	206	! Allocate types for all but last iteration
	207	call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
	208	gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels,Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
	209	gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
	210	gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
	211	gbx%use_precipitation_fluxes,gbx%use_reff, &
	212	gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
	213	gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
	214	gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
	215	gbx_it)
	216	call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
	217	call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
	218	call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
	219	call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
	220	call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
	221	call construct_cosp_misr(cfg,Ni,misr_it)
	222	call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
	223	call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
	224	elseif (i == Niter) then ! last iteration
	225	call free_cosp_gridbox(gbx_it,.true.)
	226	call free_cosp_subgrid(sgx_it)
	227	call free_cosp_vgrid(vgrid_it)
	228	call free_cosp_sgradar(sgradar_it)
	229	call free_cosp_sglidar(sglidar_it)
	230	call free_cosp_isccp(isccp_it)
	231	call free_cosp_misr(misr_it)
	232	call free_cosp_radarstats(stradar_it)
	233	call free_cosp_lidarstats(stlidar_it)
	234	! Allocate types for iterations
	235	call construct_cosp_gridbox(gbx%time,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables,gbx%use_gas_abs, &
	236	gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels,Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
	237	gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
	238	gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
	239	gbx%use_precipitation_fluxes,gbx%use_reff, &
	240	gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
	241	gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
	242	gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
	243	gbx_it)
	244	! --- Copy arrays without Npoints as dimension ---
	245	gbx_it%dist_prmts_hydro = gbx%dist_prmts_hydro
	246	gbx_it%dist_type_aero = gbx_it%dist_type_aero
	247	call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
	248	call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
	249	call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
	250	call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
	251	call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
	252	call construct_cosp_misr(cfg,Ni,misr_it)
	253	call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
	254	call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
	255	endif
	256	! --- Copy sections of arrays with Npoints as dimension ---
	257	ix=(/i_first,i_last/)
	258	iy=(/1,Ni/)
	259	call cosp_gridbox_cpsection(ix,iy,gbx,gbx_it)
	260	! These serve as initialisation of *_it types
	261	call cosp_subgrid_cpsection(ix,iy,sgx,sgx_it)
	262	if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar,sgradar_it)
	263	if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar,sglidar_it)
	264	if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp,isccp_it)
	265	if (cfg%Lmisr_sim) call cosp_misr_cpsection(ix,iy,misr,misr_it)
	266	if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar,stradar_it)
	267	if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar,stlidar_it)
	268	print *,'---------ix: ',ix
	269	call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
	270	sglidar_it,isccp_it,misr_it,stradar_it,stlidar_it)
	271
	272	! --- Copy results to output structures ---
	273	! call cosp_gridbox_cphp(gbx_it,gbx)
	274	ix=(/1,Ni/)
	275	iy=(/i_first,i_last/)
	276	call cosp_subgrid_cpsection(ix,iy,sgx_it,sgx)
	277	if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar_it,sgradar)
	278	if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar_it,sglidar)
	279	if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp_it,isccp)
	280	if (cfg%Lmisr_sim) call cosp_misr_cpsection(ix,iy,misr_it,misr)
	281	if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar_it,stradar)
	282	if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar_it,stlidar)
	283	enddo
	284	! Deallocate types
	285	call free_cosp_gridbox(gbx_it,.true.)
	286	call free_cosp_subgrid(sgx_it)
	287	call free_cosp_vgrid(vgrid_it)
	288	call free_cosp_sgradar(sgradar_it)
	289	call free_cosp_sglidar(sglidar_it)
	290	call free_cosp_isccp(isccp_it)
	291	call free_cosp_misr(misr_it)
	292	call free_cosp_radarstats(stradar_it)
	293	call free_cosp_lidarstats(stlidar_it)
	294	endif
	295
	296
	297	END SUBROUTINE COSP
	298
	299	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	300	!--------------------- SUBROUTINE COSP_ITER ----------------------
	301	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	302	SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
	303
	304	! Arguments
	305	integer,intent(in) :: overlap ! overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
	306	integer,dimension(:),intent(in) :: seed
	307	type(cosp_config),intent(in) :: cfg ! Configuration options
	308	type(cosp_vgrid),intent(in) :: vgrid ! Information on vertical grid of stats
	309	type(cosp_gridbox),intent(inout) :: gbx
	310	type(cosp_subgrid),intent(inout) :: sgx ! Subgrid info
	311	type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
	312	type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
	313	type(cosp_isccp),intent(inout) :: isccp ! Output from ISCCP simulator
	314	type(cosp_misr),intent(inout) :: misr ! Output from MISR simulator
	315	type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
	316	type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
	317
	318	! Local variables
	319	integer :: Npoints ! Number of gridpoints
	320	integer :: Ncolumns ! Number of subcolumns
	321	integer :: Nlevels ! Number of levels
	322	integer :: Nhydro ! Number of hydrometeors
	323	integer :: Niter ! Number of calls to cosp_simulator
	324	integer :: i,j,k
[1414]	325	integer :: I_HYDRO
[1262]	326	real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
[1327]	327	integer,parameter :: scops_debug=0 ! set to non-zero value to print out inputs for debugging in SCOPS
	328
[1262]	329	real,dimension(:, :),allocatable :: cca_scops,ls_p_rate,cv_p_rate, &
	330	tca_scops ! Cloud cover in each model level (HORIZONTAL gridbox fraction) of total cloud.
	331	! Levels are from TOA to SURFACE. (nPoints, nLev)
	332	real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
	333	! Levels are from SURFACE to TOA
[1414]	334	real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric dens
[1262]	335	type(cosp_sghydro) :: sghydro ! Subgrid info for hydrometeors en each iteration
	336
	337
	338	!++++++++++ Dimensions ++++++++++++
	339	Npoints = gbx%Npoints
	340	Ncolumns = gbx%Ncolumns
	341	Nlevels = gbx%Nlevels
	342	Nhydro = gbx%Nhydro
	343
	344
	345	!++++++++++ Climate/NWP mode ++++++++++
	346	if (Ncolumns > 1) then
	347	!++++++++++ Subgrid sampling ++++++++++
	348	! Allocate arrays before calling SCOPS
	349	allocate(frac_ls(Npoints,Nlevels),frac_cv(Npoints,Nlevels),prec_ls(Npoints,Nlevels),prec_cv(Npoints,Nlevels))
	350	allocate(tca_scops(Npoints,Nlevels),cca_scops(Npoints,Nlevels), &
	351	ls_p_rate(Npoints,Nlevels),cv_p_rate(Npoints,Nlevels))
	352	! Initialize to zero
	353	frac_ls=0.0
	354	prec_ls=0.0
	355	frac_cv=0.0
	356	prec_cv=0.0
	357	! Cloud fractions for SCOPS from TOA to SFC
	358	tca_scops = gbx%tca(:,Nlevels:1:-1)
	359	cca_scops = gbx%cca(:,Nlevels:1:-1)
	360
	361	! Call to SCOPS
	362	! strat and conv arrays are passed with levels from TOA to SURFACE.
	363	call scops(Npoints,Nlevels,Ncolumns,seed,tca_scops,cca_scops,overlap,sgx%frac_out,scops_debug)
	364
	365	! temporarily use prec_ls/cv to transfer information about precipitation flux into prec_scops
	366	if(gbx%use_precipitation_fluxes) then
	367	ls_p_rate(:,Nlevels:1:-1)=gbx%rain_ls(:,1:Nlevels)+gbx%snow_ls(:,1:Nlevels)+gbx%grpl_ls(:,1:Nlevels)
	368	cv_p_rate(:,Nlevels:1:-1)=gbx%rain_cv(:,1:Nlevels)+gbx%snow_cv(:,1:Nlevels)
	369	else
	370	ls_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_LSRAIN)+ &
	371	gbx%mr_hydro(:,1:Nlevels,I_LSSNOW)+ &
	372	gbx%mr_hydro(:,1:Nlevels,I_LSGRPL)
	373	cv_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_CVRAIN)+ &
	374	gbx%mr_hydro(:,1:Nlevels,I_CVSNOW)
	375	endif
	376
	377	call prec_scops(Npoints,Nlevels,Ncolumns,ls_p_rate,cv_p_rate,sgx%frac_out,sgx%prec_frac)
	378
	379	! Precipitation fraction
	380	do j=1,Npoints,1
	381	do k=1,Nlevels,1
	382	do i=1,Ncolumns,1
[1414]	383	if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
	384	if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
[1262]	385	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
	386	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
	387	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 3) then
	388	prec_cv(j,k)=prec_cv(j,k)+1.
	389	prec_ls(j,k)=prec_ls(j,k)+1.
	390	endif
	391	enddo !i
	392	frac_ls(j,k)=frac_ls(j,k)/Ncolumns
	393	frac_cv(j,k)=frac_cv(j,k)/Ncolumns
	394	prec_ls(j,k)=prec_ls(j,k)/Ncolumns
	395	prec_cv(j,k)=prec_cv(j,k)/Ncolumns
	396	enddo !k
	397	enddo !j
	398
	399	! Levels from SURFACE to TOA.
	400	if (NpointsNcolumnsNlevels < 10000) then
	401	sgx%frac_out(1:Npoints,:,1:Nlevels) = sgx%frac_out(1:Npoints,:,Nlevels:1:-1)
	402	sgx%prec_frac(1:Npoints,:,1:Nlevels) = sgx%prec_frac(1:Npoints,:,Nlevels:1:-1)
	403	else
	404	! This is done within a loop (unvectorized) over nPoints to save memory
	405	do j=1,Npoints
	406	sgx%frac_out(j,:,1:Nlevels) = sgx%frac_out(j,:,Nlevels:1:-1)
	407	sgx%prec_frac(j,:,1:Nlevels) = sgx%prec_frac(j,:,Nlevels:1:-1)
	408	enddo
	409	endif
	410
	411	! Deallocate arrays that will no longer be used
	412	deallocate(tca_scops,cca_scops,ls_p_rate,cv_p_rate)
	413
	414	! Populate the subgrid arrays
	415	call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
	416	do k=1,Ncolumns
	417	!--------- Mixing ratios for clouds and Reff for Clouds and precip -------
	418	column_frac_out => sgx%frac_out(:,k,:)
[1414]	419	where (column_frac_out == I_LSC) !+++++++++++ LS clouds ++++++++
[1262]	420	sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
	421	sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
	422
	423	sghydro%Reff(:,k,:,I_LSCLIQ) = gbx%Reff(:,:,I_LSCLIQ)
	424	sghydro%Reff(:,k,:,I_LSCICE) = gbx%Reff(:,:,I_LSCICE)
	425	sghydro%Reff(:,k,:,I_LSRAIN) = gbx%Reff(:,:,I_LSRAIN)
	426	sghydro%Reff(:,k,:,I_LSSNOW) = gbx%Reff(:,:,I_LSSNOW)
	427	sghydro%Reff(:,k,:,I_LSGRPL) = gbx%Reff(:,:,I_LSGRPL)
[1414]	428	elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++
[1262]	429	sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ)
	430	sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE)
	431
	432	sghydro%Reff(:,k,:,I_CVCLIQ) = gbx%Reff(:,:,I_CVCLIQ)
	433	sghydro%Reff(:,k,:,I_CVCICE) = gbx%Reff(:,:,I_CVCICE)
	434	sghydro%Reff(:,k,:,I_CVRAIN) = gbx%Reff(:,:,I_CVRAIN)
	435	sghydro%Reff(:,k,:,I_CVSNOW) = gbx%Reff(:,:,I_CVSNOW)
	436	end where
	437	!--------- Precip -------
	438	if (.not. gbx%use_precipitation_fluxes) then
[1414]	439	where (column_frac_out == I_LSC) !+++++++++++ LS Precipitation ++++++++
[1262]	440	sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN)
	441	sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW)
	442	sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL)
[1414]	443	elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++
[1262]	444	sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN)
	445	sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW)
	446	end where
	447	endif
	448	enddo
	449	! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values
	450	do k=1,Nlevels
	451	do j=1,Npoints
	452	!--------- Clouds -------
	453	if (frac_ls(j,k) .ne. 0.) then
	454	sghydro%mr_hydro(j,:,k,I_LSCLIQ) = sghydro%mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
	455	sghydro%mr_hydro(j,:,k,I_LSCICE) = sghydro%mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
	456	endif
	457	if (frac_cv(j,k) .ne. 0.) then
	458	sghydro%mr_hydro(j,:,k,I_CVCLIQ) = sghydro%mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
	459	sghydro%mr_hydro(j,:,k,I_CVCICE) = sghydro%mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
	460	endif
	461	!--------- Precip -------
	462	if (gbx%use_precipitation_fluxes) then
	463	if (prec_ls(j,k) .ne. 0.) then
	464	gbx%rain_ls(j,k) = gbx%rain_ls(j,k)/prec_ls(j,k)
	465	gbx%snow_ls(j,k) = gbx%snow_ls(j,k)/prec_ls(j,k)
	466	gbx%grpl_ls(j,k) = gbx%grpl_ls(j,k)/prec_ls(j,k)
	467	endif
	468	if (prec_cv(j,k) .ne. 0.) then
	469	gbx%rain_cv(j,k) = gbx%rain_cv(j,k)/prec_cv(j,k)
	470	gbx%snow_cv(j,k) = gbx%snow_cv(j,k)/prec_cv(j,k)
	471	endif
	472	else
	473	if (prec_ls(j,k) .ne. 0.) then
	474	sghydro%mr_hydro(j,:,k,I_LSRAIN) = sghydro%mr_hydro(j,:,k,I_LSRAIN)/prec_ls(j,k)
	475	sghydro%mr_hydro(j,:,k,I_LSSNOW) = sghydro%mr_hydro(j,:,k,I_LSSNOW)/prec_ls(j,k)
	476	sghydro%mr_hydro(j,:,k,I_LSGRPL) = sghydro%mr_hydro(j,:,k,I_LSGRPL)/prec_ls(j,k)
	477	endif
	478	if (prec_cv(j,k) .ne. 0.) then
	479	sghydro%mr_hydro(j,:,k,I_CVRAIN) = sghydro%mr_hydro(j,:,k,I_CVRAIN)/prec_cv(j,k)
	480	sghydro%mr_hydro(j,:,k,I_CVSNOW) = sghydro%mr_hydro(j,:,k,I_CVSNOW)/prec_cv(j,k)
	481	endif
	482	endif
	483	enddo !k
	484	enddo !j
	485	deallocate(frac_ls,prec_ls,frac_cv,prec_cv)
	486
	487	if (gbx%use_precipitation_fluxes) then
	488	! convert precipitation flux into mixing ratio
	489	call pf_to_mr(Npoints,Nlevels,Ncolumns,gbx%rain_ls,gbx%snow_ls,gbx%grpl_ls, &
	490	gbx%rain_cv,gbx%snow_cv,sgx%prec_frac,gbx%p,gbx%T, &
	491	sghydro%mr_hydro(:,:,:,I_LSRAIN),sghydro%mr_hydro(:,:,:,I_LSSNOW),sghydro%mr_hydro(:,:,:,I_LSGRPL), &
	492	sghydro%mr_hydro(:,:,:,I_CVRAIN),sghydro%mr_hydro(:,:,:,I_CVSNOW))
[1414]	493	endif
[1262]	494	!++++++++++ CRM mode ++++++++++
	495	else
	496	sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
	497	sghydro%Reff(:,1,:,:) = gbx%Reff
	498	!--------- Clouds -------
	499	where ((gbx%dtau_s > 0.0))
	500	sgx%frac_out(:,1,:) = 1 ! Subgrid cloud array. Dimensions (Npoints,Ncolumns,Nlevels)
	501	endwhere
	502	endif ! Ncolumns > 1
	503
	504
	505	!++++++++++ Simulator ++++++++++
	506	call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,stradar,stlidar)
	507
	508	! Deallocate subgrid arrays
	509	call free_cosp_sghydro(sghydro)
	510	END SUBROUTINE COSP_ITER
	511
	512	END MODULE MOD_COSP

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

Download in other formats:

Original Format