Context Navigation

scops.F @ 5453

Last change on this file since 5453 was 1262, checked in by idelkadi, 15 years ago
-Rajout des fichiers cosp_input_nl.txt et cosp_output_nl.txt pour controler les parametres entrees sorties COSP -Rajout du repertoire cosp contenant les sources de COSP
File size: 11.4 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	implicit none
42
43	INTEGER npoints ! number of model points in the horizontal
44	INTEGER nlev ! number of model levels in column
45	INTEGER ncol ! number of subcolumns
46
47
48	INTEGER overlap ! overlap type
49	! 1=max
50	! 2=rand
51	! 3=max/rand
52	REAL cc(npoints,nlev)
53	! input cloud cover in each model level (fraction)
54	! NOTE: This is the HORIZONTAL area of each
55	! grid box covered by clouds
56
57	REAL conv(npoints,nlev)
58	! input convective cloud cover in each model
59	! level (fraction)
60	! NOTE: This is the HORIZONTAL area of each
61	! grid box covered by convective clouds
62
63	INTEGER i,j,ilev,ibox,ncolprint,ilev2
64
65	REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
66	! Equivalent of BOX in original version, but
67	! indexed by column then row, rather than
68	! by row then column
69
70
71	INTEGER seed(npoints)
72	! seed values for marsaglia random number generator
73	! It is recommended that the seed is set
74	! to a different value for each model
75	! gridbox it is called on, as it is
76	! possible that the choice of the same
77	! seed value every time may introduce some
78	! statistical bias in the results, particularly
79	! for low values of NCOL.
80
81	REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction)
82	! with extra layer of zeroes on top
83	! in this version this just contains the values input
84	! from cc but with an extra level
85
86	REAL threshold(npoints,ncol) ! pointer to position in gridbox
87	REAL maxocc(npoints,ncol) ! Flag for max overlapped conv cld
88	REAL maxosc(npoints,ncol) ! Flag for max overlapped strat cld
89
90	REAL boxpos(npoints,ncol) ! ordered pointer to position in gridbox
91
92	REAL threshold_min(npoints,ncol) ! minimum value to define range in with new threshold
93	! is chosen
94
95	REAL ran(npoints) ! vector of random numbers
96
97	INTEGER irand,i2_16,huge32,overflow_32 ! variables for RNG
98	PARAMETER(huge32=2147483647)
99	i2_16=65536
100
101	do ibox=1,ncol
102	do j=1,npoints
103	boxpos(j,ibox)=(ibox-.5)/ncol
104	enddo
105	enddo
106
107	! ---------------------------------------------------!
108	! Initialise working variables
109	! ---------------------------------------------------!
110
111	! Initialised frac_out to zero
112
113	do ilev=1,nlev
114	do ibox=1,ncol
115	do j=1,npoints
116	frac_out(j,ibox,ilev)=0.0
117	enddo
118	enddo
119	enddo
120
121	! assign 2d tca array using 1d input array cc
122
123	do j=1,npoints
124	tca(j,0)=0
125	enddo
126
127	do ilev=1,nlev
128	do j=1,npoints
129	tca(j,ilev)=cc(j,ilev)
130	enddo
131	enddo
132
133	if (ncolprint.ne.0) then
134	write (6,'(a)') 'frac_out_pp_rev:'
135	do j=1,npoints,1000
136	write(6,'(a10)') 'j='
137	write(6,'(8I10)') j
138	write (6,'(8f5.2)')
139	& ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
140
141	enddo
142	write (6,'(a)') 'ncol:'
143	write (6,'(I3)') ncol
144	endif
145	if (ncolprint.ne.0) then
146	write (6,'(a)') 'last_frac_pp:'
147	do j=1,npoints,1000
148	write(6,'(a10)') 'j='
149	write(6,'(8I10)') j
150	write (6,'(8f5.2)') (tca(j,0))
151	enddo
152	endif
153
154	! ---------------------------------------------------!
155	! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
156	! frac_out is the array that contains the information
157	! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
158	! convective cloud
159
160	!loop over vertical levels
161	DO 200 ilev = 1,nlev
162
163	! Initialise threshold
164
165	IF (ilev.eq.1) then
166	! If max overlap
167	IF (overlap.eq.1) then
168	! select pixels spread evenly
169	! across the gridbox
170	DO ibox=1,ncol
171	do j=1,npoints
172	threshold(j,ibox)=boxpos(j,ibox)
173	enddo
174	enddo
175	ELSE
176	DO ibox=1,ncol
177	include 'congvec.h'
178	! select random pixels from the non-convective
179	! part the gridbox ( some will be converted into
180	! convective pixels below )
181	do j=1,npoints
182	threshold(j,ibox)=
183	& conv(j,ilev)+(1-conv(j,ilev))*ran(j)
184	enddo
185	enddo
186	ENDIF
187	IF (ncolprint.ne.0) then
188	write (6,'(a)') 'threshold_nsf2:'
189	do j=1,npoints,1000
190	write(6,'(a10)') 'j='
191	write(6,'(8I10)') j
192	write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
193	enddo
194	ENDIF
195	ENDIF
196
197	IF (ncolprint.ne.0) then
198	write (6,'(a)') 'ilev:'
199	write (6,'(I2)') ilev
200	ENDIF
201
202	DO ibox=1,ncol
203
204	! All versions
205	do j=1,npoints
206	if (boxpos(j,ibox).le.conv(j,ilev)) then
207	maxocc(j,ibox) = 1.
208	else
209	maxocc(j,ibox) = 0.
210	end if
211	enddo
212
213	! Max overlap
214	if (overlap.eq.1) then
215	do j=1,npoints
216	threshold_min(j,ibox)=conv(j,ilev)
217	maxosc(j,ibox)=1
218	enddo
219	endif
220
221	! Random overlap
222	if (overlap.eq.2) then
223	do j=1,npoints
224	threshold_min(j,ibox)=conv(j,ilev)
225	maxosc(j,ibox)=0
226	enddo
227	endif
228
229	! Max/Random overlap
230	if (overlap.eq.3) then
231	do j=1,npoints
232	threshold_min(j,ibox)=max(conv(j,ilev),
233	& min(tca(j,ilev-1),tca(j,ilev)))
234	if (threshold(j,ibox)
235	& .lt.min(tca(j,ilev-1),tca(j,ilev))
236	& .and.(threshold(j,ibox).gt.conv(j,ilev))) then
237	maxosc(j,ibox)= 1
238	else
239	maxosc(j,ibox)= 0
240	end if
241	enddo
242	endif
243
244	! Reset threshold
245
246	include 'congvec.h'
247
248	do j=1,npoints
249	threshold(j,ibox)=
250	!if max overlapped conv cloud
251	& maxocc(j,ibox) * (
252	& boxpos(j,ibox)
253	& ) +
254	!else
255	& (1-maxocc(j,ibox)) * (
256	!if max overlapped strat cloud
257	& (maxosc(j,ibox)) * (
258	!threshold=boxpos
259	& threshold(j,ibox)
260	& ) +
261	!else
262	& (1-maxosc(j,ibox)) * (
263	!threshold_min=random[thrmin,1]
264	& threshold_min(j,ibox)+
265	& (1-threshold_min(j,ibox))*ran(j)
266	& )
267	& )
268	enddo
269
270	ENDDO ! ibox
271
272	! Fill frac_out with 1's where tca is greater than the threshold
273
274	DO ibox=1,ncol
275	do j=1,npoints
276	if (tca(j,ilev).gt.threshold(j,ibox)) then
277	frac_out(j,ibox,ilev)=1
278	else
279	frac_out(j,ibox,ilev)=0
280	end if
281	enddo
282	ENDDO
283
284	! Code to partition boxes into startiform and convective parts
285	! goes here
286
287	DO ibox=1,ncol
288	do j=1,npoints
289	if (threshold(j,ibox).le.conv(j,ilev)) then
290	! = 2 IF threshold le conv(j)
291	frac_out(j,ibox,ilev) = 2
292	else
293	! = the same IF NOT threshold le conv(j)
294	frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev)
295	end if
296	enddo
297	ENDDO
298
299	! Set last_frac to tca at this level, so as to be tca
300	! from last level next time round
301
302	if (ncolprint.ne.0) then
303
304	do j=1,npoints ,1000
305	write(6,'(a10)') 'j='
306	write(6,'(8I10)') j
307	write (6,'(a)') 'last_frac:'
308	write (6,'(8f5.2)') (tca(j,ilev-1))
309
310	write (6,'(a)') 'conv:'
311	write (6,'(8f5.2)') (conv(j,ilev),ibox=1,ncolprint)
312
313	write (6,'(a)') 'max_overlap_cc:'
314	write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
315
316	write (6,'(a)') 'max_overlap_sc:'
317	write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
318
319	write (6,'(a)') 'threshold_min_nsf2:'
320	write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
321
322	write (6,'(a)') 'threshold_nsf2:'
323	write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
324
325	write (6,'(a)') 'frac_out_pp_rev:'
326	write (6,'(8f5.2)')
327	& ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
328	enddo
329	endif
330
331	200 CONTINUE !loop over nlev
332
333
334	end
335

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

Download in other formats:

Original Format