source: LMDZ5/trunk/libf/phylmd/cosp/scops.F @ 2884

Last change on this file since 2884 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp
  • 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
333200   CONTINUE    !loop over nlev
334
335
336      end
337
Note: See TracBrowser for help on using the repository browser.