Context Navigation

mod_cosp.F90 @ 5408

Last change on this file since 5408 was 5271, checked in by abarral, 8 weeks ago
Move dimensions.h into a module Nb: doesn't compile yet
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 31.3 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
[5271]	25	INCLUDE "cosp_defs.h"
[1262]	26	MODULE MOD_COSP
	27	USE MOD_COSP_TYPES
	28	USE MOD_COSP_SIMULATOR
[2428]	29	USE MOD_COSP_MODIS_SIMULATOR
[1262]	30	IMPLICIT NONE
	31
	32	CONTAINS
	33
	34
	35	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	36	!--------------------- SUBROUTINE COSP ---------------------------
	37	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]	38	!#ifdef RTTOV
	39	!SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
	40	!#else
	41	SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
	42	!#endif
[1262]	43	! Arguments
	44	integer,intent(in) :: overlap ! overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
	45	integer,intent(in) :: Ncolumns
	46	type(cosp_config),intent(in) :: cfg ! Configuration options
	47	type(cosp_vgrid),intent(in) :: vgrid ! Information on vertical grid of stats
	48	type(cosp_gridbox),intent(inout) :: gbx
	49	type(cosp_subgrid),intent(inout) :: sgx ! Subgrid info
	50	type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
	51	type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
	52	type(cosp_isccp),intent(inout) :: isccp ! Output from ISCCP simulator
	53	type(cosp_misr),intent(inout) :: misr ! Output from MISR simulator
[2428]	54	type(cosp_modis),intent(inout) :: modis ! Output from MODIS simulator
	55	!#ifdef RTTOV
	56	! type(cosp_rttov),intent(inout) :: rttov ! Output from RTTOV
	57	!#endif
[1262]	58	type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
	59	type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
	60
	61	! Local variables
	62	integer :: Npoints ! Number of gridpoints
	63	integer :: Nlevels ! Number of levels
	64	integer :: Nhydro ! Number of hydrometeors
	65	integer :: Niter ! Number of calls to cosp_simulator
	66	integer :: i_first,i_last ! First and last gridbox to be processed in each iteration
[2428]	67	integer :: i,Ni
[1262]	68	integer,dimension(2) :: ix,iy
	69	logical :: reff_zero
	70	real :: maxp,minp
[2428]	71	integer,dimension(:),allocatable :: & ! Dimensions nPoints
[1262]	72	seed ! It is recommended that the seed is set to a different value for each model
	73	! gridbox it is called on, as it is possible that the choice of the same
	74	! seed value every time may introduce some statistical bias in the results,
	75	! particularly for low values of NCOL.
	76	! Types used in one iteration
	77	type(cosp_gridbox) :: gbx_it
	78	type(cosp_subgrid) :: sgx_it
	79	type(cosp_vgrid) :: vgrid_it
	80	type(cosp_sgradar) :: sgradar_it
	81	type(cosp_sglidar) :: sglidar_it
	82	type(cosp_isccp) :: isccp_it
[2428]	83	type(cosp_modis) :: modis_it
[1262]	84	type(cosp_misr) :: misr_it
[2428]	85	!#ifdef RTTOV
	86	! type(cosp_rttov) :: rttov_it
	87	!#endif
[1262]	88	type(cosp_radarstats) :: stradar_it
	89	type(cosp_lidarstats) :: stlidar_it
[2428]	90
	91	!++++++++++ Dimensions ++++++++++++
[1262]	92	Npoints = gbx%Npoints
	93	Nlevels = gbx%Nlevels
	94	Nhydro = gbx%Nhydro
	95
[2428]	96	!++++++++++ Depth of model layers ++++++++++++
	97	do i=1,Nlevels-1
	98	gbx%dlev(:,i) = gbx%zlev_half(:,i+1) - gbx%zlev_half(:,i)
	99	enddo
	100	gbx%dlev(:,Nlevels) = 2.0*(gbx%zlev(:,Nlevels) - gbx%zlev_half(:,Nlevels))
	101
[1262]	102	!++++++++++ Apply sanity checks to inputs ++++++++++
	103	! call cosp_check_input('longitude',gbx%longitude,min_val=0.0,max_val=360.0)
	104	call cosp_check_input('longitude',gbx%longitude,min_val=-180.0,max_val=180.0)
	105	call cosp_check_input('latitude',gbx%latitude,min_val=-90.0,max_val=90.0)
	106	call cosp_check_input('dlev',gbx%dlev,min_val=0.0)
	107	call cosp_check_input('p',gbx%p,min_val=0.0)
	108	call cosp_check_input('ph',gbx%ph,min_val=0.0)
	109	call cosp_check_input('T',gbx%T,min_val=0.0)
	110	call cosp_check_input('q',gbx%q,min_val=0.0)
	111	call cosp_check_input('sh',gbx%sh,min_val=0.0)
	112	call cosp_check_input('dtau_s',gbx%dtau_s,min_val=0.0)
	113	call cosp_check_input('dtau_c',gbx%dtau_c,min_val=0.0)
	114	call cosp_check_input('dem_s',gbx%dem_s,min_val=0.0,max_val=1.0)
	115	call cosp_check_input('dem_c',gbx%dem_c,min_val=0.0,max_val=1.0)
	116	! Point information (Npoints)
	117	call cosp_check_input('land',gbx%land,min_val=0.0,max_val=1.0)
	118	call cosp_check_input('psfc',gbx%psfc,min_val=0.0)
	119	call cosp_check_input('sunlit',gbx%sunlit,min_val=0.0,max_val=1.0)
	120	call cosp_check_input('skt',gbx%skt,min_val=0.0)
	121	! TOTAL and CONV cloud fraction for SCOPS
	122	call cosp_check_input('tca',gbx%tca,min_val=0.0,max_val=1.0)
	123	call cosp_check_input('cca',gbx%cca,min_val=0.0,max_val=1.0)
	124	! Precipitation fluxes on model levels
	125	call cosp_check_input('rain_ls',gbx%rain_ls,min_val=0.0)
	126	call cosp_check_input('rain_cv',gbx%rain_cv,min_val=0.0)
	127	call cosp_check_input('snow_ls',gbx%snow_ls,min_val=0.0)
	128	call cosp_check_input('snow_cv',gbx%snow_cv,min_val=0.0)
	129	call cosp_check_input('grpl_ls',gbx%grpl_ls,min_val=0.0)
	130	! Hydrometeors concentration and distribution parameters
	131	call cosp_check_input('mr_hydro',gbx%mr_hydro,min_val=0.0)
	132	! Effective radius [m]. (Npoints,Nlevels,Nhydro)
	133	call cosp_check_input('Reff',gbx%Reff,min_val=0.0)
	134	reff_zero=.true.
	135	if (any(gbx%Reff > 1.e-8)) then
	136	reff_zero=.false.
	137	! reff_zero == .false.
	138	! and gbx%use_reff == .true. Reff use in radar and lidar
	139	! and reff_zero == .false. Reff use in lidar and set to 0 for radar
	140	endif
	141	if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
	142	gbx%Reff = DEFAULT_LIDAR_REFF
	143	print *, '---------- COSP WARNING ------------'
	144	print *, ''
	145	print *, 'Using default Reff in lidar simulations'
	146	print *, ''
	147	print *, '----------------------------------'
	148	endif
	149
	150	! Aerosols concentration and distribution parameters
	151	call cosp_check_input('conc_aero',gbx%conc_aero,min_val=0.0)
	152	! Checks for CRM mode
	153	if (Ncolumns == 1) then
	154	if (gbx%use_precipitation_fluxes) then
	155	print *, '---------- COSP ERROR ------------'
	156	print *, ''
	157	print *, 'Use of precipitation fluxes not supported in CRM mode (Ncolumns=1)'
	158	print *, ''
	159	print *, '----------------------------------'
	160	stop
	161	endif
	162	if ((maxval(gbx%dtau_c) > 0.0).or.(maxval(gbx%dem_c) > 0.0)) then
	163	print *, '---------- COSP ERROR ------------'
	164	print *, ''
	165	print *, ' dtau_c > 0.0 or dem_c > 0.0. In CRM mode (Ncolumns=1), '
	166	print *, ' the optical depth (emmisivity) of all clouds must be '
	167	print *, ' passed through dtau_s (dem_s)'
	168	print *, ''
	169	print *, '----------------------------------'
	170	stop
	171	endif
	172	endif
	173
	174	! We base the seed in the decimal part of the surface pressure.
[2428]	175	allocate(seed(Npoints))
	176	seed = int(gbx%psfc) ! This is to avoid division by zero when Npoints = 1
[1262]	177	! Roj Oct/2008 ... Note: seed value of 0 caused me some problems + I want to
	178	! randomize for each call to COSP even when Npoints ==1
[2428]	179	minp = minval(gbx%psfc)
	180	maxp = maxval(gbx%psfc)
	181	if (Npoints .gt. 1) seed=int((gbx%psfc-minp)/(maxp-minp)*100000) + 1
	182	! Below it's how it was done in the original implementation of the ISCCP simulator.
	183	! The one above is better for offline data, when you may have packed data
	184	! that subsamples the decimal fraction of the surface pressure.
	185	! if (Npoints .gt. 1) seed=(gbx%psfc-int(gbx%psfc))*1000000
[1262]	186
[2428]	187
[1262]	188	if (gbx%Npoints_it >= gbx%Npoints) then ! One iteration gbx%Npoints
[2428]	189	!#ifdef RTTOV
	190	! call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
	191	!#else
	192	call cosp_iter(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
	193	!#endif
[1262]	194	else ! Several iterations to save memory
	195	Niter = gbx%Npoints/gbx%Npoints_it ! Integer division
	196	if (Niter*gbx%Npoints_it < gbx%Npoints) Niter = Niter + 1
	197	do i=1,Niter
	198	i_first = (i-1)*gbx%Npoints_it + 1
	199	i_last = i_first + gbx%Npoints_it - 1
	200	i_last = min(i_last,gbx%Npoints)
	201	Ni = i_last - i_first + 1
	202	if (i == 1) then
	203	! Allocate types for all but last iteration
[2428]	204	call construct_cosp_gridbox(gbx%time,gbx%time_bnds,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables, &
	205	gbx%use_gas_abs,gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels, &
	206	Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
[1262]	207	gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
	208	gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
	209	gbx%use_precipitation_fluxes,gbx%use_reff, &
	210	gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
	211	gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
	212	gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
	213	gbx_it)
	214	call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
	215	call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
	216	call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
	217	call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
	218	call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
[2428]	219	call construct_cosp_modis(cfg, Ni, modis_it)
[1262]	220	call construct_cosp_misr(cfg,Ni,misr_it)
[2428]	221	!#ifdef RTTOV
	222	! call construct_cosp_rttov(Ni,gbx%nchan,rttov_it)
	223	!#endif
[1262]	224	call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
	225	call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
	226	elseif (i == Niter) then ! last iteration
	227	call free_cosp_gridbox(gbx_it,.true.)
	228	call free_cosp_subgrid(sgx_it)
	229	call free_cosp_vgrid(vgrid_it)
	230	call free_cosp_sgradar(sgradar_it)
	231	call free_cosp_sglidar(sglidar_it)
	232	call free_cosp_isccp(isccp_it)
[2428]	233	call free_cosp_modis(modis_it)
[1262]	234	call free_cosp_misr(misr_it)
[2428]	235	!#ifdef RTTOV
	236	! call free_cosp_rttov(rttov_it)
	237	!#endif
[1262]	238	call free_cosp_radarstats(stradar_it)
	239	call free_cosp_lidarstats(stlidar_it)
	240	! Allocate types for iterations
[2428]	241	call construct_cosp_gridbox(gbx%time,gbx%time_bnds,gbx%radar_freq,gbx%surface_radar,gbx%use_mie_tables, &
	242	gbx%use_gas_abs,gbx%do_ray,gbx%melt_lay,gbx%k2,Ni,Nlevels, &
	243	Ncolumns,N_HYDRO,gbx%Nprmts_max_hydro, &
[1262]	244	gbx%Naero,gbx%Nprmts_max_aero,Ni,gbx%lidar_ice_type,gbx%isccp_top_height, &
	245	gbx%isccp_top_height_direction,gbx%isccp_overlap,gbx%isccp_emsfc_lw, &
	246	gbx%use_precipitation_fluxes,gbx%use_reff, &
	247	gbx%plat,gbx%sat,gbx%inst,gbx%nchan,gbx%ZenAng, &
	248	gbx%Ichan(1:gbx%nchan),gbx%surfem(1:gbx%nchan), &
	249	gbx%co2,gbx%ch4,gbx%n2o,gbx%co, &
	250	gbx_it)
	251	! --- Copy arrays without Npoints as dimension ---
	252	gbx_it%dist_prmts_hydro = gbx%dist_prmts_hydro
	253	gbx_it%dist_type_aero = gbx_it%dist_type_aero
	254	call construct_cosp_vgrid(gbx_it,vgrid%Nlvgrid,vgrid%use_vgrid,vgrid%csat_vgrid,vgrid_it)
	255	call construct_cosp_subgrid(Ni, Ncolumns, Nlevels, sgx_it)
	256	call construct_cosp_sgradar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,sgradar_it)
	257	call construct_cosp_sglidar(cfg,Ni,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_it)
	258	call construct_cosp_isccp(cfg,Ni,Ncolumns,Nlevels,isccp_it)
[2428]	259	call construct_cosp_modis(cfg,Ni, modis_it)
[1262]	260	call construct_cosp_misr(cfg,Ni,misr_it)
[2428]	261	!#ifdef RTTOV
	262	! call construct_cosp_rttov(Ni,gbx%nchan,rttov_it)
	263	!#endif
[1262]	264	call construct_cosp_radarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_it)
	265	call construct_cosp_lidarstats(cfg,Ni,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_it)
	266	endif
	267	! --- Copy sections of arrays with Npoints as dimension ---
	268	ix=(/i_first,i_last/)
	269	iy=(/1,Ni/)
	270	call cosp_gridbox_cpsection(ix,iy,gbx,gbx_it)
	271	! These serve as initialisation of *_it types
	272	call cosp_subgrid_cpsection(ix,iy,sgx,sgx_it)
	273	if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar,sgradar_it)
	274	if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar,sglidar_it)
	275	if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp,isccp_it)
[2428]	276	if (cfg%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis,modis_it)
[1262]	277	if (cfg%Lmisr_sim) call cosp_misr_cpsection(ix,iy,misr,misr_it)
[2428]	278	!#ifdef RTTOV
	279	! if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov,rttov_it)
	280	!#endif
[1262]	281	if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar,stradar_it)
	282	if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar,stlidar_it)
[2428]	283	!#ifdef RTTOV
	284	! call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
	285	! sglidar_it,isccp_it,misr_it,modis_it,rttov_it,stradar_it,stlidar_it)
	286	!#else
[1262]	287	call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
[2428]	288	sglidar_it,isccp_it,misr_it,modis_it,stradar_it,stlidar_it)
	289	!#endif
[1262]	290	! --- Copy results to output structures ---
	291	ix=(/1,Ni/)
	292	iy=(/i_first,i_last/)
	293	call cosp_subgrid_cpsection(ix,iy,sgx_it,sgx)
	294	if (cfg%Lradar_sim) call cosp_sgradar_cpsection(ix,iy,sgradar_it,sgradar)
	295	if (cfg%Llidar_sim) call cosp_sglidar_cpsection(ix,iy,sglidar_it,sglidar)
	296	if (cfg%Lisccp_sim) call cosp_isccp_cpsection(ix,iy,isccp_it,isccp)
[2428]	297	if (cfg%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis_it,modis)
[1262]	298	if (cfg%Lmisr_sim) call cosp_misr_cpsection(ix,iy,misr_it,misr)
[2428]	299	!#ifdef RTTOV
	300	! if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov_it,rttov)
	301	!#endif
[1262]	302	if (cfg%Lradar_sim) call cosp_radarstats_cpsection(ix,iy,stradar_it,stradar)
	303	if (cfg%Llidar_sim) call cosp_lidarstats_cpsection(ix,iy,stlidar_it,stlidar)
	304	enddo
	305	! Deallocate types
	306	call free_cosp_gridbox(gbx_it,.true.)
	307	call free_cosp_subgrid(sgx_it)
	308	call free_cosp_vgrid(vgrid_it)
	309	call free_cosp_sgradar(sgradar_it)
	310	call free_cosp_sglidar(sglidar_it)
	311	call free_cosp_isccp(isccp_it)
[2428]	312	call free_cosp_modis(modis_it)
[1262]	313	call free_cosp_misr(misr_it)
[2428]	314	!#ifdef RTTOV
	315	! call free_cosp_rttov(rttov_it)
	316	!#endif
[1262]	317	call free_cosp_radarstats(stradar_it)
	318	call free_cosp_lidarstats(stlidar_it)
	319	endif
[2428]	320	deallocate(seed)
[1262]	321
	322
	323	END SUBROUTINE COSP
	324
	325	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	326	!--------------------- SUBROUTINE COSP_ITER ----------------------
	327	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]	328	!#ifdef RTTOV
	329	!SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
	330	!#else
	331	SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
	332	!#endif
[1262]	333	! Arguments
	334	integer,intent(in) :: overlap ! overlap type in SCOPS: 1=max, 2=rand, 3=max/rand
	335	integer,dimension(:),intent(in) :: seed
	336	type(cosp_config),intent(in) :: cfg ! Configuration options
	337	type(cosp_vgrid),intent(in) :: vgrid ! Information on vertical grid of stats
	338	type(cosp_gridbox),intent(inout) :: gbx
	339	type(cosp_subgrid),intent(inout) :: sgx ! Subgrid info
	340	type(cosp_sgradar),intent(inout) :: sgradar ! Output from radar simulator
	341	type(cosp_sglidar),intent(inout) :: sglidar ! Output from lidar simulator
	342	type(cosp_isccp),intent(inout) :: isccp ! Output from ISCCP simulator
	343	type(cosp_misr),intent(inout) :: misr ! Output from MISR simulator
[2428]	344	type(cosp_modis),intent(inout) :: modis ! Output from MODIS simulator
	345	!#ifdef RTTOV
	346	! type(cosp_rttov),intent(inout) :: rttov ! Output from RTTOV
	347	!#endif
[1262]	348	type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics from radar simulator
	349	type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics from lidar simulator
	350
	351	! Local variables
	352	integer :: Npoints ! Number of gridpoints
	353	integer :: Ncolumns ! Number of subcolumns
	354	integer :: Nlevels ! Number of levels
	355	integer :: Nhydro ! Number of hydrometeors
	356	integer :: i,j,k
[2428]	357	integer :: I_HYDRO
[1262]	358	real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
[2428]	359	real,dimension(:,:),pointer :: column_prec_out ! Array with one column of prec_frac
	360	integer :: scops_debug=0 ! set to non-zero value to print out inputs for debugging in SCOPS
[1262]	361	real,dimension(:, :),allocatable :: cca_scops,ls_p_rate,cv_p_rate, &
	362	tca_scops ! Cloud cover in each model level (HORIZONTAL gridbox fraction) of total cloud.
	363	! Levels are from TOA to SURFACE. (nPoints, nLev)
	364	real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
	365	! Levels are from SURFACE to TOA
[2428]	366	real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric density
[1262]	367	type(cosp_sghydro) :: sghydro ! Subgrid info for hydrometeors en each iteration
	368
	369
	370	!++++++++++ Dimensions ++++++++++++
	371	Npoints = gbx%Npoints
	372	Ncolumns = gbx%Ncolumns
	373	Nlevels = gbx%Nlevels
	374	Nhydro = gbx%Nhydro
	375
	376	!++++++++++ Climate/NWP mode ++++++++++
	377	if (Ncolumns > 1) then
	378	!++++++++++ Subgrid sampling ++++++++++
	379	! Allocate arrays before calling SCOPS
	380	allocate(frac_ls(Npoints,Nlevels),frac_cv(Npoints,Nlevels),prec_ls(Npoints,Nlevels),prec_cv(Npoints,Nlevels))
	381	allocate(tca_scops(Npoints,Nlevels),cca_scops(Npoints,Nlevels), &
	382	ls_p_rate(Npoints,Nlevels),cv_p_rate(Npoints,Nlevels))
	383	! Initialize to zero
	384	frac_ls=0.0
	385	prec_ls=0.0
	386	frac_cv=0.0
	387	prec_cv=0.0
	388	! Cloud fractions for SCOPS from TOA to SFC
	389	tca_scops = gbx%tca(:,Nlevels:1:-1)
	390	cca_scops = gbx%cca(:,Nlevels:1:-1)
	391
	392	! Call to SCOPS
	393	! strat and conv arrays are passed with levels from TOA to SURFACE.
	394	call scops(Npoints,Nlevels,Ncolumns,seed,tca_scops,cca_scops,overlap,sgx%frac_out,scops_debug)
	395
	396	! temporarily use prec_ls/cv to transfer information about precipitation flux into prec_scops
	397	if(gbx%use_precipitation_fluxes) then
	398	ls_p_rate(:,Nlevels:1:-1)=gbx%rain_ls(:,1:Nlevels)+gbx%snow_ls(:,1:Nlevels)+gbx%grpl_ls(:,1:Nlevels)
	399	cv_p_rate(:,Nlevels:1:-1)=gbx%rain_cv(:,1:Nlevels)+gbx%snow_cv(:,1:Nlevels)
	400	else
	401	ls_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_LSRAIN)+ &
	402	gbx%mr_hydro(:,1:Nlevels,I_LSSNOW)+ &
	403	gbx%mr_hydro(:,1:Nlevels,I_LSGRPL)
	404	cv_p_rate(:,Nlevels:1:-1)=gbx%mr_hydro(:,1:Nlevels,I_CVRAIN)+ &
	405	gbx%mr_hydro(:,1:Nlevels,I_CVSNOW)
	406	endif
	407
	408	call prec_scops(Npoints,Nlevels,Ncolumns,ls_p_rate,cv_p_rate,sgx%frac_out,sgx%prec_frac)
	409
	410	! Precipitation fraction
	411	do j=1,Npoints,1
	412	do k=1,Nlevels,1
	413	do i=1,Ncolumns,1
[1414]	414	if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
	415	if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
[1262]	416	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
	417	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
	418	if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 3) then
	419	prec_cv(j,k)=prec_cv(j,k)+1.
	420	prec_ls(j,k)=prec_ls(j,k)+1.
	421	endif
	422	enddo !i
	423	frac_ls(j,k)=frac_ls(j,k)/Ncolumns
	424	frac_cv(j,k)=frac_cv(j,k)/Ncolumns
	425	prec_ls(j,k)=prec_ls(j,k)/Ncolumns
	426	prec_cv(j,k)=prec_cv(j,k)/Ncolumns
	427	enddo !k
	428	enddo !j
	429
	430	! Levels from SURFACE to TOA.
	431	if (NpointsNcolumnsNlevels < 10000) then
	432	sgx%frac_out(1:Npoints,:,1:Nlevels) = sgx%frac_out(1:Npoints,:,Nlevels:1:-1)
	433	sgx%prec_frac(1:Npoints,:,1:Nlevels) = sgx%prec_frac(1:Npoints,:,Nlevels:1:-1)
	434	else
	435	! This is done within a loop (unvectorized) over nPoints to save memory
	436	do j=1,Npoints
	437	sgx%frac_out(j,:,1:Nlevels) = sgx%frac_out(j,:,Nlevels:1:-1)
	438	sgx%prec_frac(j,:,1:Nlevels) = sgx%prec_frac(j,:,Nlevels:1:-1)
	439	enddo
	440	endif
	441
	442	! Deallocate arrays that will no longer be used
	443	deallocate(tca_scops,cca_scops,ls_p_rate,cv_p_rate)
[2428]	444
[1262]	445	! Populate the subgrid arrays
	446	call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
	447	do k=1,Ncolumns
	448	!--------- Mixing ratios for clouds and Reff for Clouds and precip -------
	449	column_frac_out => sgx%frac_out(:,k,:)
[1414]	450	where (column_frac_out == I_LSC) !+++++++++++ LS clouds ++++++++
[1262]	451	sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
	452	sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
[2428]	453
[1262]	454	sghydro%Reff(:,k,:,I_LSCLIQ) = gbx%Reff(:,:,I_LSCLIQ)
	455	sghydro%Reff(:,k,:,I_LSCICE) = gbx%Reff(:,:,I_LSCICE)
[2428]	456
	457	sghydro%Np(:,k,:,I_LSCLIQ) = gbx%Np(:,:,I_LSCLIQ)
	458	sghydro%Np(:,k,:,I_LSCICE) = gbx%Np(:,:,I_LSCICE)
	459
[1414]	460	elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++
[2428]	461	sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ)
	462	sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE)
	463
	464	sghydro%Reff(:,k,:,I_CVCLIQ) = gbx%Reff(:,:,I_CVCLIQ)
	465	sghydro%Reff(:,k,:,I_CVCICE) = gbx%Reff(:,:,I_CVCICE)
	466
	467	sghydro%Np(:,k,:,I_CVCLIQ) = gbx%Np(:,:,I_CVCLIQ)
	468	sghydro%Np(:,k,:,I_CVCICE) = gbx%Np(:,:,I_CVCICE)
	469
	470	end where
	471	column_prec_out => sgx%prec_frac(:,k,:)
	472	where ((column_prec_out == 1) .or. (column_prec_out == 3) ) !++++ LS precip ++++
	473	sghydro%Reff(:,k,:,I_LSRAIN) = gbx%Reff(:,:,I_LSRAIN)
	474	sghydro%Reff(:,k,:,I_LSSNOW) = gbx%Reff(:,:,I_LSSNOW)
	475	sghydro%Reff(:,k,:,I_LSGRPL) = gbx%Reff(:,:,I_LSGRPL)
	476
	477	sghydro%Np(:,k,:,I_LSRAIN) = gbx%Np(:,:,I_LSRAIN)
	478	sghydro%Np(:,k,:,I_LSSNOW) = gbx%Np(:,:,I_LSSNOW)
	479	sghydro%Np(:,k,:,I_LSGRPL) = gbx%Np(:,:,I_LSGRPL)
	480	elsewhere ((column_prec_out == 2) .or. (column_prec_out == 3)) !++++ CONV precip ++++
	481	sghydro%Reff(:,k,:,I_CVRAIN) = gbx%Reff(:,:,I_CVRAIN)
	482	sghydro%Reff(:,k,:,I_CVSNOW) = gbx%Reff(:,:,I_CVSNOW)
	483
	484	sghydro%Np(:,k,:,I_CVRAIN) = gbx%Np(:,:,I_CVRAIN)
	485	sghydro%Np(:,k,:,I_CVSNOW) = gbx%Np(:,:,I_CVSNOW)
	486	end where
[1262]	487	!--------- Precip -------
	488	if (.not. gbx%use_precipitation_fluxes) then
[1414]	489	where (column_frac_out == I_LSC) !+++++++++++ LS Precipitation ++++++++
[1262]	490	sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN)
	491	sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW)
	492	sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL)
[1414]	493	elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++
[2428]	494	sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN)
	495	sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW)
[1262]	496	end where
	497	endif
	498	enddo
	499	! convert the mixing ratio and precipitation flux from gridbox mean to the fraction-based values
	500	do k=1,Nlevels
	501	do j=1,Npoints
	502	!--------- Clouds -------
	503	if (frac_ls(j,k) .ne. 0.) then
	504	sghydro%mr_hydro(j,:,k,I_LSCLIQ) = sghydro%mr_hydro(j,:,k,I_LSCLIQ)/frac_ls(j,k)
	505	sghydro%mr_hydro(j,:,k,I_LSCICE) = sghydro%mr_hydro(j,:,k,I_LSCICE)/frac_ls(j,k)
	506	endif
	507	if (frac_cv(j,k) .ne. 0.) then
	508	sghydro%mr_hydro(j,:,k,I_CVCLIQ) = sghydro%mr_hydro(j,:,k,I_CVCLIQ)/frac_cv(j,k)
	509	sghydro%mr_hydro(j,:,k,I_CVCICE) = sghydro%mr_hydro(j,:,k,I_CVCICE)/frac_cv(j,k)
	510	endif
	511	!--------- Precip -------
	512	if (gbx%use_precipitation_fluxes) then
	513	if (prec_ls(j,k) .ne. 0.) then
	514	gbx%rain_ls(j,k) = gbx%rain_ls(j,k)/prec_ls(j,k)
	515	gbx%snow_ls(j,k) = gbx%snow_ls(j,k)/prec_ls(j,k)
	516	gbx%grpl_ls(j,k) = gbx%grpl_ls(j,k)/prec_ls(j,k)
	517	endif
	518	if (prec_cv(j,k) .ne. 0.) then
	519	gbx%rain_cv(j,k) = gbx%rain_cv(j,k)/prec_cv(j,k)
	520	gbx%snow_cv(j,k) = gbx%snow_cv(j,k)/prec_cv(j,k)
	521	endif
	522	else
	523	if (prec_ls(j,k) .ne. 0.) then
	524	sghydro%mr_hydro(j,:,k,I_LSRAIN) = sghydro%mr_hydro(j,:,k,I_LSRAIN)/prec_ls(j,k)
	525	sghydro%mr_hydro(j,:,k,I_LSSNOW) = sghydro%mr_hydro(j,:,k,I_LSSNOW)/prec_ls(j,k)
	526	sghydro%mr_hydro(j,:,k,I_LSGRPL) = sghydro%mr_hydro(j,:,k,I_LSGRPL)/prec_ls(j,k)
	527	endif
	528	if (prec_cv(j,k) .ne. 0.) then
	529	sghydro%mr_hydro(j,:,k,I_CVRAIN) = sghydro%mr_hydro(j,:,k,I_CVRAIN)/prec_cv(j,k)
	530	sghydro%mr_hydro(j,:,k,I_CVSNOW) = sghydro%mr_hydro(j,:,k,I_CVSNOW)/prec_cv(j,k)
	531	endif
	532	endif
	533	enddo !k
	534	enddo !j
	535	deallocate(frac_ls,prec_ls,frac_cv,prec_cv)
	536
	537	if (gbx%use_precipitation_fluxes) then
[2428]	538
	539	#ifdef MMF_V3p5_TWO_MOMENT
	540
	541	write(,) 'Precipitation Flux to Mixing Ratio conversion not (yet?) supported ', &
	542	'for MMF3.5 Two Moment Microphysics'
	543	stop
	544	#else
	545	! Density
	546	allocate(rho(Npoints,Nlevels))
	547	I_HYDRO = I_LSRAIN
	548	call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
	549	n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
	550	g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
	551	gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
	552	gbx%rain_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
	553	I_HYDRO = I_LSSNOW
	554	call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
	555	n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
	556	g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
	557	gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
	558	gbx%snow_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
	559	I_HYDRO = I_CVRAIN
	560	call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,2., &
	561	n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
	562	g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
	563	gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
	564	gbx%rain_cv,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
	565	I_HYDRO = I_CVSNOW
	566	call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,2., &
	567	n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
	568	g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
	569	gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
	570	gbx%snow_cv,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
	571	I_HYDRO = I_LSGRPL
	572	call cosp_precip_mxratio(Npoints,Nlevels,Ncolumns,gbx%p,gbx%T,sgx%prec_frac,1., &
	573	n_ax(I_HYDRO),n_bx(I_HYDRO),alpha_x(I_HYDRO),c_x(I_HYDRO),d_x(I_HYDRO), &
	574	g_x(I_HYDRO),a_x(I_HYDRO),b_x(I_HYDRO), &
	575	gamma_1(I_HYDRO),gamma_2(I_HYDRO),gamma_3(I_HYDRO),gamma_4(I_HYDRO), &
	576	gbx%grpl_ls,sghydro%mr_hydro(:,:,:,I_HYDRO),sghydro%Reff(:,:,:,I_HYDRO))
	577	if(allocated(rho)) deallocate(rho)
	578	#endif
	579
	580	endif
[1262]	581	!++++++++++ CRM mode ++++++++++
	582	else
[2428]	583	call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
[1262]	584	sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
	585	sghydro%Reff(:,1,:,:) = gbx%Reff
[2428]	586	sghydro%Np(:,1,:,:) = gbx%Np ! added by Roj with Quickbeam V3.0
	587
[1262]	588	!--------- Clouds -------
	589	where ((gbx%dtau_s > 0.0))
	590	sgx%frac_out(:,1,:) = 1 ! Subgrid cloud array. Dimensions (Npoints,Ncolumns,Nlevels)
	591	endwhere
	592	endif ! Ncolumns > 1
	593
	594	!++++++++++ Simulator ++++++++++
[2428]	595	!#ifdef RTTOV
	596	! call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
	597	!#else
	598	call cosp_simulator(gbx,sgx,sghydro,cfg,vgrid,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
	599	!#endif
[1262]	600
	601	! Deallocate subgrid arrays
	602	call free_cosp_sghydro(sghydro)
	603	END SUBROUTINE COSP_ITER
	604
	605	END MODULE MOD_COSP

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

Download in other formats:

Original Format