source: LMDZ4/branches/LMDZ4-dev/libf/cosp/scops.F @ 1269

Last change on this file since 1269 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
331200   CONTINUE    !loop over nlev
332
333
334      end
335
Note: See TracBrowser for help on using the repository browser.