source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/scrip/src/fracnnei.f @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 4 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818

File size: 9.4 KB
Line 
1      subroutine fracnnei (src_size, dst_size,
2     $    ld_srcmask, ld_dstmask,
3     $    src_lon, src_lat, dst_lon, dst_lat,
4     $    num_links, num_wgts,
5     $    weights, src_addr, dst_addr)
6C****
7C               *****************************
8C               * OASIS ROUTINE  -  LEVEL 4 *
9C               * -------------     ------- *
10C               *****************************
11C
12C**** *fracnnei* - SCRIP remapping
13C
14C     Purpose:
15C     -------
16C     Treatment of the tricky points in an interpolation
17C
18C     Interface:
19C     ---------
20C       *CALL*  *
21C
22C     Called from:
23C     -----------
24C     scriprmp
25C
26C     Input:
27C     -----
28C             src_size    : source grid size (integer)
29C             dst_size    : target grid size (integer)
30C             ld_srcmask  : mask of the source grid
31C             ld_dstmask  : mask of the target grid
32C             src_lon     : longitudes of the source grid
33C             src_lat     : latitudes of the source grid
34C             dst_lon     : longitudes of the target grid
35C             dst_lat     : latitudes of the target grid
36C             num_links   : total number of links
37C             num_wgts    : number of weights for each link
38C     InOut
39C     -----
40C             weights     : remapping weights
41C             src_addr    : remapping source addresses
42C             dst_addr    : remapping target addresses
43C
44C     History:
45C     -------
46C       Version   Programmer     Date        Description
47C       -------   ----------     ----        ----------- 
48C       2.5       D.Declat       2002/08/20  adapted from S. Valcke ptmsq
49C       3.0       S. Valcke      2002/10/30  test and corrections
50C
51C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52C* ---------------------------- Modules used ----------------------------
53C
54      use kinds_mod     ! defines common data types
55      use constants     ! defines common constants
56      use grids         ! module containing grid information
57      use remap_vars    ! module containing remap information
58      USE mod_oasis_flush
59C
60C* ---------------------------- Implicit --------------------------------
61C
62      implicit none
63C
64C* ---------------------------- Include files ---------------------------
65C
66
67C      INCLUDE 'netcdf.inc'
68C
69C* ---------------------------- Intent In -------------------------------
70C
71      INTEGER (kind=int_kind) ::
72     $    src_size,             ! size of the source grid
73     $    dst_size              ! size of the destination grid
74C
75      REAL (kind=dbl_kind) ::
76     $    src_lat(src_size), src_lon(src_size),
77     $    dst_lat(dst_size), dst_lon(dst_size)
78
79C
80      LOGICAL ::
81     $    ld_srcmask(src_size),   ! source grid mask
82     $    ld_dstmask(dst_size)    ! target grid mask
83C
84      INTEGER (kind=int_kind) ::
85     $    num_links,      ! number of links between src and tgt
86     $    num_wgts        ! number of weights
87
88C
89C* ---------------------------- Intent InOut ------------------------------
90C
91      REAL (kind=dbl_kind) ::
92     $    weights(num_wgts, num_links) ! remapping weights
93C
94      INTEGER (kind=int_kind) ::
95     $    src_addr(num_links), ! remapping source addresses
96     $    dst_addr(num_links)  ! remapping target addresses
97C
98C* ---------------------------- Local declarations ----------------------
99C
100C
101      INTEGER (kind=int_kind) :: 
102     $    ila_nneiadd           ! Nearest-neighbor address
103C
104      INTEGER (kind=int_kind) ::
105     $    ib_dst,               ! INDEX loop for the distance grid
106     $    ib_src,               ! INDEX loop for the source grid
107     $    ib_neigh,             ! INDEX loop on the corresponding src pts
108     $    ib_links              ! INDEX loop for the links     
109C
110      INTEGER (kind=int_kind) ::
111     $    nb_Vmm,               ! number of Mtt points
112     $    ntottotland,          ! number of land points
113     $    ntotland,             ! number of land points
114     $    ntotoce,              ! number of oceanic points
115     $    ntotngh               ! number of corresponding source points
116C
117      INTEGER (kind=int_kind) ::
118     $    beg_links,            ! begining of the serie of links
119     $    nb_links              ! number of links
120C
121      REAL (kind=dbl_kind) ::
122     $    coslat,               ! cosinus of the latitude
123     $    sinlat,               ! sinus of the latitude
124     $    coslon,               ! cosinus of the longitude
125     $    sinlon,               ! sinus of the longitude
126     $    distance, 
127     $    dist_min,
128     $    arg
129C
130      INTEGER (kind=int_kind) :: n, il
131C
132C* ---------------------------- Poema verses ----------------------------
133C
134C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135C
136C*    1. Initialization
137C        --------------
138C
139      IF (nlogprt .GE. 2) THEN
140          WRITE (UNIT = nulou,FMT = *) ' '
141          WRITE (UNIT = nulou,FMT = *) ' '
142          WRITE (UNIT = nulou,FMT = *) 
143     $        '   Entering ROUTINE fracnnei  -  Level 4'
144          WRITE (UNIT = nulou,FMT = *) 
145     $        '           ****************     *******'
146          WRITE (UNIT = nulou,FMT = *) ' '
147          WRITE (UNIT = nulou,FMT = *) 
148     $        ' Treating the tricky points of the remapping'
149          WRITE (UNIT = nulou,FMT = *) ' '
150          CALL OASIS_FLUSH_SCRIP(nulou)
151      ENDIF
152C
153C *----------------------------------------------------------------------
154C
155C*    2. Treating Vmm points   V
156C        -------------------  m m
157C     The target point is a non-masked Valid point while the source points 
158C         are all masked points. Use of the non-masked nearest neighbours.
159C
160      nb_Vmm = 0
161      ntottotland = 0
162      beg_links = 1
163
164C* -- Loop all other the target points
165      DO ib_dst = 1, dst_size
166C* -- If the point is a sea point
167        IF (ld_dstmask(ib_dst)) THEN
168
169            beg_links = 0
170
171            DO ib_links = 1, num_links
172              IF (dst_addr(ib_links) .eq. ib_dst) THEN
173                  beg_links = ib_links
174                  exit
175              ENDIF
176            END DO
177
178            IF (beg_links .ne. 0) THEN
179
180                ntotland = 0
181                ntotoce = 0
182                ntotngh = 0
183
184C* -- Find the number of associated src points to the non-masked tgt point
185                nb_links = 1
186                DO il = beg_links+1, num_links
187                  IF (dst_addr(il) .eq. dst_addr(beg_links)) 
188     $                nb_links = nb_links + 1
189                END DO
190               
191C* -- For each point on the src grid associated to the non-masked tgt point
192                DO ib_neigh = 1, nb_links
193C
194                  ntotngh = ntotngh + 1
195
196C* -- Check IF the point is a masked or non-masked point and treat it
197                  IF (.not. ld_srcmask(src_addr(beg_links+ib_neigh-1)))
198     $                THEN
199                      ntotland = ntotland + 1
200                      ntottotland = ntottotland + 1
201                  ELSE IF (ld_srcmask(src_addr(beg_links+ib_neigh-1))) 
202     $                    THEN
203                      ntotoce = ntotoce + 1
204                  ELSE
205                      WRITE (nulou, *) 'Pb with ocean mask with Mtt 1' 
206                  END IF
207
208                END DO
209 
210C* -- If all the src points are land, treat it !
211                IF (ntotland .EQ. ntotngh) THEN
212                    nb_Vmm = nb_Vmm + 1
213                    WRITE(nulou,*) 
214     $               '************ Doing FRACNNEI for point', ib_dst 
215
216C* -- Find the nearest neighbours and change weights and address 
217
218                    coslat = cos(dst_lat(ib_dst))
219                    sinlat = sin(dst_lat(ib_dst))
220                    coslon = cos(dst_lon(ib_dst))
221                    sinlon = sin(dst_lon(ib_dst))
222
223                    dist_min = bignum
224                    ila_nneiadd = 0
225                    DO ib_src = 1, src_size
226                      IF (ld_srcmask(ib_src)) THEN
227                          arg = 
228     &                        coslat*cos(src_lat(ib_src))*
229     &                       (coslon*cos(src_lon(ib_src)) +
230     &                        sinlon*sin(src_lon(ib_src)))+
231     &                        sinlat*sin(src_lat(ib_src))
232                          IF (arg < -1.0d0) THEN
233                              arg = -1.0d0
234                          ELSE IF (arg > 1.0d0) THEN
235                              arg = 1.0d0
236                          END IF
237                          distance = acos(arg)
238                          IF (distance < dist_min) THEN
239                              ila_nneiadd = ib_src
240                              dist_min = distance
241                          ENDIF
242                      ENDIF
243                    END DO
244                    src_addr(beg_links) = ila_nneiadd
245                    weights(1,beg_links) = 1.0
246                    WRITE(nulou,*) 
247     $               '*************** Nearest source neighbour is ', 
248     $                  ila_nneiadd                   
249                    IF (nb_links > 1) then
250                        DO n=2, nb_links
251                          src_addr(beg_links+n-1)=0
252                          weights(1,beg_links+n-1)=0.0
253                        END DO
254                    END IF
255                END IF
256            END IF
257        ENDIF
258      END DO
259C
260C
261C *----------------------------------------------------------------------
262C
263      IF (nlogprt .GE. 2) THEN
264          WRITE (UNIT = nulou,FMT = *) ' '
265          WRITE (UNIT = nulou,FMT = *) 
266     $        '   Leaving ROUTINE fracnnei  -  Level 4'
267          WRITE (UNIT = nulou,FMT = *) ' '
268          CALL OASIS_FLUSH_SCRIP(nulou)
269      ENDIF
270
271      END SUBROUTINE fracnnei
272
273!***********************************************************************
274
275 
Note: See TracBrowser for help on using the repository browser.