Context Navigation

scops.F @ 5409

Last change on this file since 5409 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

Rev	Line
[1262]	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.
[2428]	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 $
[5099]	10
[1262]	11	! Redistribution and use in source and binary forms, with or without
	12	! modification, are permitted provided that the
	13	! following conditions are met:
[5099]	14
[1262]	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.
[5099]	26
[1262]	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.
[5099]	38
[1262]	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
[5095]	135	if (ncolprint.ne.0) then
[1262]	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
[5095]	147	if (ncolprint.ne.0) then
[1262]	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
[5095]	163	DO 200 ilev = 1,nlev
[1262]	164
	165	! Initialise threshold
	166
[5095]	167	IF (ilev.eq.1) then
[1262]	168	! If max overlap
[5095]	169	IF (overlap.eq.1) then
[1262]	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
[2428]	179	include 'congvec.h'
[1262]	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
[5095]	189	IF (ncolprint.ne.0) then
[1262]	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
[5095]	199	IF (ncolprint.ne.0) then
[1262]	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
[5095]	208	if (boxpos(j,ibox).le.conv(j,ilev)) then
[2428]	209	maxocc(j,ibox) = 1
[1262]	210	else
[2428]	211	maxocc(j,ibox) = 0
[1262]	212	end if
	213	enddo
	214
	215	! Max overlap
[5095]	216	if (overlap.eq.1) then
[1262]	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
[5095]	224	if (overlap.eq.2) then
[1262]	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
[5095]	232	if (overlap.eq.3) then
[1262]	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)
[5095]	237	& .lt.min(tca(j,ilev-1),tca(j,ilev))
[5185]	238	& .AND.(threshold(j,ibox).gt.conv(j,ilev))) then
[1262]	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
[5095]	278	if (tca(j,ilev).gt.threshold(j,ibox)) then
[1262]	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
[5095]	291	if (threshold(j,ibox).le.conv(j,ilev)) then
[1262]	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
[5095]	304	if (ncolprint.ne.0) then
[1262]	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
[5095]	333	200 CONTINUE !loop over nlev
[1262]	334
	335
	336	end
	337

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

Download in other formats:

Original Format