source: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90 @ 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
File size: 10.9 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2009, British Crown Copyright, the Met Office
3! All rights reserved.
4
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28
29! History
30! May 2015 - D. Swales - Modified for COSPv2.0
31! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32module mod_scops
33  USE COSP_KINDS,     ONLY: wp
34  USE MOD_RNG!,        ONLY: rng_state,get_rng
35  use mod_cosp_error, ONLY: errorMessage
36
37  implicit none
38
39  integer,parameter :: default_overlap = 3 ! Used when invalid overlap assumption is provided.
40 
41contains
42  subroutine scops(npoints,nlev,ncol,rngs,cc,conv,overlap,frac_out,ncolprint)
43    INTEGER :: npoints,    &    ! Number of model points in the horizontal
44               nlev,       &    ! Number of model levels in column
45               ncol,       &    ! Number of subcolumns
46               overlap          ! Overlap type (1=max, 2=rand, 3=max/rand)
47    type(rng_state),dimension(npoints) :: rngs           
48    INTEGER, parameter :: huge32 = 2147483647
49    INTEGER, parameter :: i2_16  = 65536
50    INTEGER :: i,j,ilev,ibox,ncolprint,ilev2
51
52    REAL(WP), dimension(npoints,nlev) ::  &
53         cc,         &    ! Input cloud cover in each model level (fraction)
54                          ! NOTE:  This is the HORIZONTAL area of each
55                          !        grid box covered by clouds
56         conv,       &    ! Input convective cloud cover in each model level (fraction)
57                          ! NOTE:  This is the HORIZONTAL area of each
58                          !        grid box covered by convective clouds
59         tca              ! Total cloud cover in each model level (fraction)
60                          ! with extra layer of zeroes on top
61                          ! in this version this just contains the values input
62                          ! from cc but with an extra level
63    REAL(WP),intent(inout), dimension(npoints,ncol,nlev) :: &
64         frac_out         ! Boxes gridbox divided up into equivalent of BOX in
65                          ! original version, but indexed by column then row, rather than
66                          ! by row then column
67    REAL(WP), dimension(npoints,ncol) :: &
68         threshold,   &   ! pointer to position in gridbox
69         maxocc,      &   ! Flag for max overlapped conv cld
70         maxosc,      &   ! Flag for max overlapped strat cld
71         boxpos,      &   ! ordered pointer to position in gridbox
72         threshold_min    ! minimum value to define range in with new threshold is chosen.
73    REAL(WP), dimension(npoints) :: &
74         ran              ! vector of random numbers
75
76    ! Test for valid input overlap assumption
77    if (overlap .ne. 1 .AND. overlap .ne. 2 .AND. overlap .ne. 3) then
78       overlap=default_overlap
79       call errorMessage('ERROR(scops): Invalid overlap assumption provided. Using default overlap assumption (max/ran)')
80    endif
81
82    boxpos = spread(([(i, i=1,ncol)]-0.5)/ncol,1,npoints)
83   
84    ! #######################################################################
85    ! Initialize working variables
86    ! #######################################################################
87   
88    ! Initialize frac_out to zero
89    frac_out(1:npoints,1:ncol,1:nlev)=0.0     
90   
91    ! Assign 2d tca array using 1d input array cc
92    tca(1:npoints,1:nlev) = cc(1:npoints,1:nlev)
93   
94    if (ncolprint.ne.0) then
95       write (6,'(a)') 'frac_out_pp_rev:'
96       DO j=1,npoints,1000
97          write(6,'(a10)') 'j='
98          write(6,'(8I10)') j
99          write (6,'(8f5.2)') ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
100       enddo
101       write (6,'(a)') 'ncol:'
102       write (6,'(I3)') ncol
103    endif
104    if (ncolprint.ne.0) then
105       write (6,'(a)') 'last_frac_pp:'
106       DO j=1,npoints,1000
107          write(6,'(a10)') 'j='
108          write(6,'(8I10)') j
109          write (6,'(8f5.2)') (tca(j,1))
110       enddo
111    endif
112   
113    ! #######################################################################
114    ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS
115    ! frac_out is the array that contains the information
116    ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a
117    ! convective cloud
118    ! #######################################################################
119   
120    ! Loop over vertical levels
121    DO ilev = 1,nlev
122       
123       ! Initialise threshold
124       IF (ilev.eq.1) then
125          ! If max overlap
126          IF (overlap.eq.1) then
127             ! Select pixels spread evenly across the gridbox
128             threshold(1:npoints,1:ncol)=boxpos(1:npoints,1:ncol)
129          ELSE
130             DO ibox=1,ncol
131                !include 'congvec.f90'
132                ran(1:npoints) = get_rng(RNGS)
133                ! select random pixels from the non-convective
134                ! part the gridbox ( some will be converted into
135                ! convective pixels below )
136                threshold(1:npoints,ibox) = conv(1:npoints,ilev)+(1-conv(1:npoints,ilev))*ran(npoints)
137             enddo
138          ENDIF
139          IF (ncolprint.ne.0) then
140             write (6,'(a)') 'threshold_nsf2:'
141             DO j=1,npoints,1000
142                write(6,'(a10)') 'j='
143                write(6,'(8I10)') j
144                write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
145             enddo
146          ENDIF
147       ENDIF
148       
149       IF (ncolprint.ne.0) then
150          write (6,'(a)') 'ilev:'
151          write (6,'(I2)') ilev
152       ENDIF
153       
154       DO ibox=1,ncol
155          ! All versions
156          !maxocc(1:npoints,ibox) = merge(1,0,boxpos(1:npoints,ibox) .le. conv(1:npoints,ilev))
157          !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox))
158          DO j=1,npoints
159             if (boxpos(j,ibox).le.conv(j,ilev)) then
160                maxocc(j,ibox) = 1
161             else
162                maxocc(j,ibox) = 0
163             end if
164          enddo
165         
166          ! Max overlap
167          if (overlap.eq.1) then
168             threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
169             maxosc(1:npoints,ibox)        = 1               
170          endif
171         
172          ! Random overlap
173          if (overlap.eq.2) then
174             threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
175             maxosc(1:npoints,ibox)        = 0
176          endif
177          ! Max/Random overlap
178          if (overlap.eq.3) then
179             ! DS2014 START: The bounds on tca are not valid when ilev=1.
180             !threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
181             !maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
182             !     min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .AND. &
183             !     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
184             if (ilev .ne. 1) then
185                threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
186                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
187                     min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .AND. &
188                     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
189             else
190                threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(0._wp,tca(1:npoints,ilev)))
191                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
192                     min(0._wp,tca(1:npoints,ilev)) .AND. &
193                     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
194             endif
195          endif
196         
197          ! Reset threshold
198          !include 'congvec.f90'
199          ran(1:npoints) = get_rng(RNGS)
200         
201          threshold(1:npoints,ibox)= maxocc(1:npoints,ibox)*(boxpos(1:npoints,ibox)) +            &
202               (1-maxocc(1:npoints,ibox))*((maxosc(1:npoints,ibox))*(threshold(1:npoints,ibox)) + &
203               (1-maxosc(1:npoints,ibox))*(threshold_min(1:npoints,ibox)+                         &
204               (1-threshold_min(1:npoints,ibox))*ran(1:npoints)))
205         
206          ! Fill frac_out with 1's where tca is greater than the threshold
207          frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev).gt.threshold(1:npoints,ibox))
208         
209          ! Code to partition boxes into startiform and convective parts goes here
210          where(threshold(1:npoints,ibox).le.conv(1:npoints,ilev) .AND. conv(1:npoints,ilev).gt.0.) frac_out(1:npoints,ibox,ilev)=2
211       ENDDO ! ibox
212       
213       
214       ! Set last_frac to tca at this level, so as to be tca from last level next time round
215       if (ncolprint.ne.0) then
216          DO j=1,npoints ,1000
217             write(6,'(a10)') 'j='
218             write(6,'(8I10)') j
219             write (6,'(a)') 'last_frac:'
220             write (6,'(8f5.2)') (tca(j,ilev))
221             write (6,'(a)') 'conv:'
222             write (6,'(8f5.2)') (conv(j,ilev),ibox=1,ncolprint)
223             write (6,'(a)') 'max_overlap_cc:'
224             write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
225             write (6,'(a)') 'max_overlap_sc:'
226             write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
227             write (6,'(a)') 'threshold_min_nsf2:'
228             write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
229             write (6,'(a)') 'threshold_nsf2:'
230             write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
231             write (6,'(a)') 'frac_out_pp_rev:'
232             write (6,'(8f5.2)') ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
233          enddo
234       endif
235       
236    enddo ! Loop over nlev
237   
238    ! END
239  end subroutine scops
240end module mod_scops
Note: See TracBrowser for help on using the repository browser.