source: trunk/libf/cosp/scops.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 11.5 KB
RevLine 
[1]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
335200   CONTINUE    !loop over nlev
336
337
338      end
339
Note: See TracBrowser for help on using the repository browser.