Line data Source code
1 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 : ! Copyright (c) 2009, British Crown Copyright, the Met Office
3 : ! All rights reserved.
4 : !
5 : ! Redistribution and use in source and binary forms, with or without modification, are
6 : ! permitted provided that the following conditions are met:
7 : !
8 : ! 1. Redistributions of source code must retain the above copyright notice, this list of
9 : ! conditions and the following disclaimer.
10 : !
11 : ! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12 : ! of conditions and the following disclaimer in the documentation and/or other
13 : ! materials provided with the distribution.
14 : !
15 : ! 3. Neither the name of the copyright holder nor the names of its contributors may be
16 : ! used to endorse or promote products derived from this software without specific prior
17 : ! written permission.
18 : !
19 : ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21 : ! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22 : ! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 : ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24 : ! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 : ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 : ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 : ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 : !
29 : ! History
30 : ! May 2015 - D. Swales - Modified for COSPv2.0
31 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32 : module mod_scops
33 : USE COSP_KINDS, ONLY: wp
34 : USE MOD_RNG!, ONLY: rng_state,get_rng
35 : use mod_cosp_error, ONLY: errorMessage
36 :
37 : implicit none
38 :
39 : integer,parameter :: default_overlap = 3 ! Used when invalid overlap assumption is provided.
40 :
41 : contains
42 9288 : subroutine scops(npoints,nlev,ncol,rngs,cc,conv,overlap,frac_out,ncolprint)
43 : INTEGER :: npoints, & ! Number of model points in the horizontal
44 : nlev, & ! Number of model levels in column
45 : ncol, & ! Number of subcolumns
46 : overlap ! Overlap type (1=max, 2=rand, 3=max/rand)
47 : type(rng_state),dimension(npoints) :: rngs
48 : INTEGER, parameter :: huge32 = 2147483647
49 : INTEGER, parameter :: i2_16 = 65536
50 : INTEGER :: i,j,ilev,ibox,ncolprint,ilev2
51 :
52 : REAL(WP), dimension(npoints,nlev) :: &
53 : cc, & ! Input cloud cover in each model level (fraction)
54 : ! NOTE: This is the HORIZONTAL area of each
55 : ! grid box covered by clouds
56 : conv, & ! Input convective cloud cover in each model level (fraction)
57 : ! NOTE: This is the HORIZONTAL area of each
58 : ! grid box covered by convective clouds
59 18576 : tca ! Total cloud cover in each model level (fraction)
60 : ! with extra layer of zeroes on top
61 : ! in this version this just contains the values input
62 : ! from cc but with an extra level
63 : REAL(WP),intent(inout), dimension(npoints,ncol,nlev) :: &
64 : frac_out ! Boxes gridbox divided up into equivalent of BOX in
65 : ! original version, but indexed by column then row, rather than
66 : ! by row then column
67 : REAL(WP), dimension(npoints,ncol) :: &
68 18576 : threshold, & ! pointer to position in gridbox
69 18576 : maxocc, & ! Flag for max overlapped conv cld
70 18576 : maxosc, & ! Flag for max overlapped strat cld
71 18576 : boxpos, & ! ordered pointer to position in gridbox
72 18576 : threshold_min ! minimum value to define range in with new threshold is chosen.
73 : REAL(WP), dimension(npoints) :: &
74 18576 : ran ! vector of random numbers
75 :
76 : ! Test for valid input overlap assumption
77 9288 : if (overlap .ne. 1 .and. overlap .ne. 2 .and. overlap .ne. 3) then
78 0 : overlap=default_overlap
79 0 : call errorMessage('ERROR(scops): Invalid overlap assumption provided. Using default overlap assumption (max/ran)')
80 : endif
81 :
82 195048 : boxpos = spread(([(i, i=1,ncol)]-0.5)/ncol,1,npoints)
83 :
84 : ! #######################################################################
85 : ! Initialize working variables
86 : ! #######################################################################
87 :
88 : ! Initialize frac_out to zero
89 131072688 : frac_out(1:npoints,1:ncol,1:nlev)=0.0
90 :
91 : ! Assign 2d tca array using 1d input array cc
92 13045968 : tca(1:npoints,1:nlev) = cc(1:npoints,1:nlev)
93 :
94 9288 : if (ncolprint.ne.0) then
95 0 : write (6,'(a)') 'frac_out_pp_rev:'
96 0 : do j=1,npoints,1000
97 0 : write(6,'(a10)') 'j='
98 0 : write(6,'(8I10)') j
99 0 : write (6,'(8f5.2)') ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
100 : enddo
101 0 : write (6,'(a)') 'ncol:'
102 0 : write (6,'(I3)') ncol
103 : endif
104 9288 : if (ncolprint.ne.0) then
105 0 : write (6,'(a)') 'last_frac_pp:'
106 0 : do j=1,npoints,1000
107 0 : write(6,'(a10)') 'j='
108 0 : write(6,'(8I10)') j
109 0 : write (6,'(8f5.2)') (tca(j,1))
110 : enddo
111 : endif
112 :
113 : ! #######################################################################
114 : ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
115 : ! frac_out is the array that contains the information
116 : ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
117 : ! convective cloud
118 : ! #######################################################################
119 :
120 : ! Loop over vertical levels
121 789480 : DO ilev = 1,nlev
122 :
123 : ! Initialise threshold
124 780192 : IF (ilev.eq.1) then
125 : ! If max overlap
126 9288 : IF (overlap.eq.1) then
127 : ! Select pixels spread evenly across the gridbox
128 0 : threshold(1:npoints,1:ncol)=boxpos(1:npoints,1:ncol)
129 : ELSE
130 102168 : DO ibox=1,ncol
131 : !include 'congvec.f90'
132 92880 : ran(1:npoints) = get_rng(RNGS)
133 : ! select random pixels from the non-convective
134 : ! part the gridbox ( some will be converted into
135 : ! convective pixels below )
136 1560168 : threshold(1:npoints,ibox) = conv(1:npoints,ilev)+(1-conv(1:npoints,ilev))*ran(npoints)
137 : enddo
138 : ENDIF
139 9288 : IF (ncolprint.ne.0) then
140 0 : write (6,'(a)') 'threshold_nsf2:'
141 0 : do j=1,npoints,1000
142 0 : write(6,'(a10)') 'j='
143 0 : write(6,'(8I10)') j
144 0 : write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
145 : enddo
146 : ENDIF
147 : ENDIF
148 :
149 780192 : IF (ncolprint.ne.0) then
150 0 : write (6,'(a)') 'ilev:'
151 0 : write (6,'(I2)') ilev
152 : ENDIF
153 :
154 8582112 : DO ibox=1,ncol
155 : ! All versions
156 : !maxocc(1:npoints,ibox) = merge(1,0,boxpos(1:npoints,ibox) .le. conv(1:npoints,ilev))
157 : !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox))
158 130273920 : do j=1,npoints
159 130273920 : if (boxpos(j,ibox).le.conv(j,ilev)) then
160 598222 : maxocc(j,ibox) = 1
161 : else
162 121873778 : maxocc(j,ibox) = 0
163 : end if
164 : enddo
165 :
166 : ! Max overlap
167 7801920 : if (overlap.eq.1) then
168 0 : threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
169 0 : maxosc(1:npoints,ibox) = 1
170 : endif
171 :
172 : ! Random overlap
173 7801920 : if (overlap.eq.2) then
174 0 : threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
175 0 : maxosc(1:npoints,ibox) = 0
176 : endif
177 : ! Max/Random overlap
178 7801920 : if (overlap.eq.3) then
179 : ! DS2014 START: The bounds on tca are not valid when ilev=1.
180 : !threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
181 : !maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
182 : ! min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. &
183 : ! (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
184 7801920 : if (ilev .ne. 1) then
185 128723040 : threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
186 7709040 : maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
187 : min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. &
188 136432080 : (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
189 : else
190 1550880 : threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(0._wp,tca(1:npoints,ilev)))
191 92880 : maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
192 : min(0._wp,tca(1:npoints,ilev)) .and. &
193 1643760 : (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
194 : endif
195 : endif
196 :
197 : ! Reset threshold
198 : !include 'congvec.f90'
199 7801920 : ran(1:npoints) = get_rng(RNGS)
200 :
201 15603840 : threshold(1:npoints,ibox)= maxocc(1:npoints,ibox)*(boxpos(1:npoints,ibox)) + &
202 : (1-maxocc(1:npoints,ibox))*((maxosc(1:npoints,ibox))*(threshold(1:npoints,ibox)) + &
203 : (1-maxosc(1:npoints,ibox))*(threshold_min(1:npoints,ibox)+ &
204 145877760 : (1-threshold_min(1:npoints,ibox))*ran(1:npoints)))
205 :
206 : ! Fill frac_out with 1's where tca is greater than the threshold
207 130273920 : frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev).gt.threshold(1:npoints,ibox))
208 :
209 : ! Code to partition boxes into startiform and convective parts goes here
210 138856032 : where(threshold(1:npoints,ibox).le.conv(1:npoints,ilev) .and. conv(1:npoints,ilev).gt.0.) frac_out(1:npoints,ibox,ilev)=2
211 : ENDDO ! ibox
212 :
213 :
214 : ! Set last_frac to tca at this level, so as to be tca from last level next time round
215 789480 : if (ncolprint.ne.0) then
216 0 : do j=1,npoints ,1000
217 0 : write(6,'(a10)') 'j='
218 0 : write(6,'(8I10)') j
219 0 : write (6,'(a)') 'last_frac:'
220 0 : write (6,'(8f5.2)') (tca(j,ilev))
221 0 : write (6,'(a)') 'conv:'
222 0 : write (6,'(8f5.2)') (conv(j,ilev),ibox=1,ncolprint)
223 0 : write (6,'(a)') 'max_overlap_cc:'
224 0 : write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
225 0 : write (6,'(a)') 'max_overlap_sc:'
226 0 : write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
227 0 : write (6,'(a)') 'threshold_min_nsf2:'
228 0 : write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
229 0 : write (6,'(a)') 'threshold_nsf2:'
230 0 : write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
231 0 : write (6,'(a)') 'frac_out_pp_rev:'
232 0 : write (6,'(8f5.2)') ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
233 : enddo
234 : endif
235 :
236 : enddo ! Loop over nlev
237 :
238 : ! END
239 9288 : end subroutine scops
240 : end module mod_scops
|
