source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/cosp/mod_cosp.F90 @ 3373

Last change on this file since 3373 was 3242, checked in by idelkadi, 7 years ago

Nettoyage du code :

  • changement de nomes des routines pour avoir les memes nome des modules
  • corrections
  • dos2unix appliquee aux fichiers
  • Rajouter "implicit none" pour les routines ne le comportant pas
  • Suppression des commentaires unitils

Corrections pour tourner en mode debug :

controle par la cles ok_debug_cosp dans mod_cosp_constants.F90
la valeur par defaut etant .false.

  • 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
Line 
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
25#include "cosp_defs.h"
26MODULE MOD_COSP
27  USE MOD_COSP_TYPES
28  USE MOD_COSP_SIMULATOR
29  USE MOD_COSP_MODIS_SIMULATOR
30  IMPLICIT NONE
31
32CONTAINS
33
34
35!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36!--------------------- SUBROUTINE COSP ---------------------------
37!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38!#ifdef RTTOV
39!SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
40!#else
41SUBROUTINE COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
42!#endif
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
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
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
67  integer :: i,Ni
68  integer,dimension(2) :: ix,iy
69  logical :: reff_zero
70  real :: maxp,minp
71  integer,dimension(:),allocatable :: & ! Dimensions nPoints
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
83  type(cosp_modis)   :: modis_it
84  type(cosp_misr)    :: misr_it
85!#ifdef RTTOV
86!  type(cosp_rttov)   :: rttov_it
87!#endif
88  type(cosp_radarstats) :: stradar_it
89  type(cosp_lidarstats) :: stlidar_it
90
91!++++++++++ Dimensions ++++++++++++
92  Npoints  = gbx%Npoints
93  Nlevels  = gbx%Nlevels
94  Nhydro   = gbx%Nhydro
95
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
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.
175   allocate(seed(Npoints))
176   seed = int(gbx%psfc) ! This is to avoid division by zero when Npoints = 1   
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
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
186
187 
188   if (gbx%Npoints_it >= gbx%Npoints) then ! One iteration gbx%Npoints
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
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
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, &
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)
219                call construct_cosp_modis(cfg, Ni, modis_it)
220                call construct_cosp_misr(cfg,Ni,misr_it)
221!#ifdef RTTOV
222!                call construct_cosp_rttov(Ni,gbx%nchan,rttov_it)
223!#endif
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)
233                call free_cosp_modis(modis_it)
234                call free_cosp_misr(misr_it)
235!#ifdef RTTOV
236!                call free_cosp_rttov(rttov_it)
237!#endif
238                call free_cosp_radarstats(stradar_it)
239                call free_cosp_lidarstats(stlidar_it)
240                ! Allocate types for iterations
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, &
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)
259                call construct_cosp_modis(cfg,Ni, modis_it)
260                call construct_cosp_misr(cfg,Ni,misr_it)
261!#ifdef RTTOV
262!                call construct_cosp_rttov(Ni,gbx%nchan,rttov_it)
263!#endif
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)
276            if (cfg%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis,modis_it)
277            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr,misr_it)
278!#ifdef RTTOV
279!            if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov,rttov_it)
280!#endif
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)
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
287            call cosp_iter(overlap,seed(ix(1):ix(2)),cfg,vgrid_it,gbx_it,sgx_it,sgradar_it, &
288                           sglidar_it,isccp_it,misr_it,modis_it,stradar_it,stlidar_it)
289!#endif
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)
297            if (cfg%Lmodis_sim) call cosp_modis_cpsection(ix,iy,modis_it,modis)
298            if (cfg%Lmisr_sim)  call cosp_misr_cpsection(ix,iy,misr_it,misr)
299!#ifdef RTTOV
300!            if (cfg%Lrttov_sim) call cosp_rttov_cpsection(ix,iy,rttov_it,rttov)
301!#endif
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)
312        call free_cosp_modis(modis_it)
313        call free_cosp_misr(misr_it)
314!#ifdef RTTOV
315!        call free_cosp_rttov(rttov_it)
316!#endif
317        call free_cosp_radarstats(stradar_it)
318        call free_cosp_lidarstats(stlidar_it)
319   endif
320   deallocate(seed)
321
322   
323END SUBROUTINE COSP
324
325!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
326!--------------------- SUBROUTINE COSP_ITER ----------------------
327!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328!#ifdef RTTOV
329!SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
330!#else
331SUBROUTINE COSP_ITER(overlap,seed,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
332!#endif
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
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
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
357  integer :: I_HYDRO
358  real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
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
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
366  real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric density
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
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.
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 (Npoints*Ncolumns*Nlevels < 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)
444
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,:)
450            where (column_frac_out == I_LSC)     !+++++++++++ LS clouds ++++++++
451                sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
452                sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
453
454                sghydro%Reff(:,k,:,I_LSCLIQ)     = gbx%Reff(:,:,I_LSCLIQ)
455                sghydro%Reff(:,k,:,I_LSCICE)     = gbx%Reff(:,:,I_LSCICE)
456
457                sghydro%Np(:,k,:,I_LSCLIQ)     = gbx%Np(:,:,I_LSCLIQ)
458                sghydro%Np(:,k,:,I_LSCICE)     = gbx%Np(:,:,I_LSCICE)
459
460            elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++
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
487            !--------- Precip -------
488            if (.not. gbx%use_precipitation_fluxes) then
489                where (column_frac_out == I_LSC)  !+++++++++++ LS Precipitation ++++++++
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)
493                elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++
494                    sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN)
495                    sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW)
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
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
581   !++++++++++ CRM mode ++++++++++
582   else
583      call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,Nhydro,sghydro)
584      sghydro%mr_hydro(:,1,:,:) = gbx%mr_hydro
585      sghydro%Reff(:,1,:,:) = gbx%Reff
586      sghydro%Np(:,1,:,:) = gbx%Np      ! added by Roj with Quickbeam V3.0
587     
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 ++++++++++
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
600
601    ! Deallocate subgrid arrays
602    call free_cosp_sghydro(sghydro)
603END SUBROUTINE COSP_ITER
604
605END MODULE MOD_COSP
Note: See TracBrowser for help on using the repository browser.