Context Navigation

cosp.F90 @ 3817

Last change on this file since 3817 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

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
[2428]	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