Context Navigation

scops.F @ 1366

Last change on this file since 1366 was 1318, checked in by yann meurdesoif, 15 years ago
YM : Parallelisation COSP MPI+OpenMP
File size: 11.5 KB

Line
1	subroutine scops(npoints,nlev,ncol,seed,cc,conv,
2	& overlap,frac_out,ncolprint)
3
4
5	! ***************************COPYRIGHT**************************
6	! (c) British Crown Copyright 2009, the Met Office.
7	! All rights reserved.
8	!
9	! Redistribution and use in source and binary forms, with or without
10	! modification, are permitted provided that the
11	! following conditions are met:
12	!
13	! * Redistributions of source code must retain the above
14	! copyright notice, this list of conditions and the following
15	! disclaimer.
16	! * Redistributions in binary form must reproduce the above
17	! copyright notice, this list of conditions and the following
18	! disclaimer in the documentation and/or other materials
19	! provided with the distribution.
20	! * Neither the name of the Met Office nor the names of its
21	! contributors may be used to endorse or promote products
22	! derived from this software without specific prior written
23	! permission.
24	!
25	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
26	! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
27	! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
28	! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
29	! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
31	! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33	! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34	! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36	!
37	! ***************************COPYRIGHT*****************************
38	! ***************************COPYRIGHT*****************************
39	! ***************************COPYRIGHT*****************************
40
41	USE mod_phys_lmdz_para
42	USE mod_grid_phy_lmdz
43
44	implicit none
45
46	INTEGER npoints ! number of model points in the horizontal
47	INTEGER nlev ! number of model levels in column
48	INTEGER ncol ! number of subcolumns
49
50
51	INTEGER overlap ! overlap type
52	! 1=max
53	! 2=rand
54	! 3=max/rand
55	REAL cc(npoints,nlev)
56	! input cloud cover in each model level (fraction)
57	! NOTE: This is the HORIZONTAL area of each
58	! grid box covered by clouds
59
60	REAL conv(npoints,nlev)
61	! input convective cloud cover in each model
62	! level (fraction)
63	! NOTE: This is the HORIZONTAL area of each
64	! grid box covered by convective clouds
65
66	INTEGER i,j,ilev,ibox,ncolprint,ilev2
67
68	REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
69	! Equivalent of BOX in original version, but
70	! indexed by column then row, rather than
71	! by row then column
72
73
74	INTEGER seed(npoints)
75	! seed values for marsaglia random number generator
76	! It is recommended that the seed is set
77	! to a different value for each model
78	! gridbox it is called on, as it is
79	! possible that the choice of the same
80	! seed value every time may introduce some
81	! statistical bias in the results, particularly
82	! for low values of NCOL.
83
84	REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction)
85	! with extra layer of zeroes on top
86	! in this version this just contains the values input
87	! from cc but with an extra level
88
89	REAL threshold(npoints,ncol) ! pointer to position in gridbox
90	REAL maxocc(npoints,ncol) ! Flag for max overlapped conv cld
91	REAL maxosc(npoints,ncol) ! Flag for max overlapped strat cld
92
93	REAL boxpos(npoints,ncol) ! ordered pointer to position in gridbox
94
95	REAL threshold_min(npoints,ncol) ! minimum value to define range in with new threshold
96	! is chosen
97
98	REAL ran(npoints) ! vector of random numbers
99
100	INTEGER irand,i2_16,huge32,overflow_32 ! variables for RNG
101	PARAMETER(huge32=2147483647)
102	i2_16=65536
103
104	do ibox=1,ncol
105	do j=1,npoints
106	boxpos(j,ibox)=(ibox-.5)/ncol
107	enddo
108	enddo
109
110	! ---------------------------------------------------!
111	! Initialise working variables
112	! ---------------------------------------------------!
113
114	! Initialised frac_out to zero
115
116	do ilev=1,nlev
117	do ibox=1,ncol
118	do j=1,npoints
119	frac_out(j,ibox,ilev)=0.0
120	enddo
121	enddo
122	enddo
123
124	! assign 2d tca array using 1d input array cc
125
126	do j=1,npoints
127	tca(j,0)=0
128	enddo
129
130	do ilev=1,nlev
131	do j=1,npoints
132	tca(j,ilev)=cc(j,ilev)
133	enddo
134	enddo
135
136	if (ncolprint.ne.0) then
137	write (6,'(a)') 'frac_out_pp_rev:'
138	do j=1,npoints,1000
139	write(6,'(a10)') 'j='
140	write(6,'(8I10)') j
141	write (6,'(8f5.2)')
142	& ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
143
144	enddo
145	write (6,'(a)') 'ncol:'
146	write (6,'(I3)') ncol
147	endif
148	if (ncolprint.ne.0) then
149	write (6,'(a)') 'last_frac_pp:'
150	do j=1,npoints,1000
151	write(6,'(a10)') 'j='
152	write(6,'(8I10)') j
153	write (6,'(8f5.2)') (tca(j,0))
154	enddo
155	endif
156
157	! ---------------------------------------------------!
158	! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
159	! frac_out is the array that contains the information
160	! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
161	! convective cloud
162
163	!loop over vertical levels
164	DO 200 ilev = 1,nlev
165
166	! Initialise threshold
167
168	IF (ilev.eq.1) then
169	! If max overlap
170	IF (overlap.eq.1) then
171	! select pixels spread evenly
172	! across the gridbox
173	DO ibox=1,ncol
174	do j=1,npoints
175	threshold(j,ibox)=boxpos(j,ibox)
176	enddo
177	enddo
178	ELSE
179	DO ibox=1,ncol
180	! include 'congvec_para.h'
181	include 'congvec.h'
182	! select random pixels from the non-convective
183	! part the gridbox ( some will be converted into
184	! convective pixels below )
185	do j=1,npoints
186	threshold(j,ibox)=
187	& conv(j,ilev)+(1-conv(j,ilev))*ran(j)
188	enddo
189	enddo
190	ENDIF
191	IF (ncolprint.ne.0) then
192	write (6,'(a)') 'threshold_nsf2:'
193	do j=1,npoints,1000
194	write(6,'(a10)') 'j='
195	write(6,'(8I10)') j
196	write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
197	enddo
198	ENDIF
199	ENDIF
200
201	IF (ncolprint.ne.0) then
202	write (6,'(a)') 'ilev:'
203	write (6,'(I2)') ilev
204	ENDIF
205
206	DO ibox=1,ncol
207
208	! All versions
209	do j=1,npoints
210	if (boxpos(j,ibox).le.conv(j,ilev)) then
211	maxocc(j,ibox) = 1.
212	else
213	maxocc(j,ibox) = 0.
214	end if
215	enddo
216
217	! Max overlap
218	if (overlap.eq.1) then
219	do j=1,npoints
220	threshold_min(j,ibox)=conv(j,ilev)
221	maxosc(j,ibox)=1
222	enddo
223	endif
224
225	! Random overlap
226	if (overlap.eq.2) then
227	do j=1,npoints
228	threshold_min(j,ibox)=conv(j,ilev)
229	maxosc(j,ibox)=0
230	enddo
231	endif
232
233	! Max/Random overlap
234	if (overlap.eq.3) then
235	do j=1,npoints
236	threshold_min(j,ibox)=max(conv(j,ilev),
237	& min(tca(j,ilev-1),tca(j,ilev)))
238	if (threshold(j,ibox)
239	& .lt.min(tca(j,ilev-1),tca(j,ilev))
240	& .and.(threshold(j,ibox).gt.conv(j,ilev))) then
241	maxosc(j,ibox)= 1
242	else
243	maxosc(j,ibox)= 0
244	end if
245	enddo
246	endif
247
248	! Reset threshold
249
250	include 'congvec.h'
251
252	do j=1,npoints
253	threshold(j,ibox)=
254	!if max overlapped conv cloud
255	& maxocc(j,ibox) * (
256	& boxpos(j,ibox)
257	& ) +
258	!else
259	& (1-maxocc(j,ibox)) * (
260	!if max overlapped strat cloud
261	& (maxosc(j,ibox)) * (
262	!threshold=boxpos
263	& threshold(j,ibox)
264	& ) +
265	!else
266	& (1-maxosc(j,ibox)) * (
267	!threshold_min=random[thrmin,1]
268	& threshold_min(j,ibox)+
269	& (1-threshold_min(j,ibox))*ran(j)
270	& )
271	& )
272	enddo
273
274	ENDDO ! ibox
275
276	! Fill frac_out with 1's where tca is greater than the threshold
277
278	DO ibox=1,ncol
279	do j=1,npoints
280	if (tca(j,ilev).gt.threshold(j,ibox)) then
281	frac_out(j,ibox,ilev)=1
282	else
283	frac_out(j,ibox,ilev)=0
284	end if
285	enddo
286	ENDDO
287
288	! Code to partition boxes into startiform and convective parts
289	! goes here
290
291	DO ibox=1,ncol
292	do j=1,npoints
293	if (threshold(j,ibox).le.conv(j,ilev)) then
294	! = 2 IF threshold le conv(j)
295	frac_out(j,ibox,ilev) = 2
296	else
297	! = the same IF NOT threshold le conv(j)
298	frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev)
299	end if
300	enddo
301	ENDDO
302
303	! Set last_frac to tca at this level, so as to be tca
304	! from last level next time round
305
306	if (ncolprint.ne.0) then
307
308	do j=1,npoints ,1000
309	write(6,'(a10)') 'j='
310	write(6,'(8I10)') j
311	write (6,'(a)') 'last_frac:'
312	write (6,'(8f5.2)') (tca(j,ilev-1))
313
314	write (6,'(a)') 'conv:'
315	write (6,'(8f5.2)') (conv(j,ilev),ibox=1,ncolprint)
316
317	write (6,'(a)') 'max_overlap_cc:'
318	write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
319
320	write (6,'(a)') 'max_overlap_sc:'
321	write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
322
323	write (6,'(a)') 'threshold_min_nsf2:'
324	write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
325
326	write (6,'(a)') 'threshold_nsf2:'
327	write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
328
329	write (6,'(a)') 'frac_out_pp_rev:'
330	write (6,'(8f5.2)')
331	& ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
332	enddo
333	endif
334
335	200 CONTINUE !loop over nlev
336
337
338	end
339

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

Download in other formats:

Original Format