Context Navigation

scops.F @ 5218

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

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

Download in other formats:

Original Format