source: TOOLS/MOSAIX/src/MOSAIX/interpol.f90 @ 4167

Last change on this file since 4167 was 4167, checked in by omamce, 4 years ago

O.M. : adapt area reading for ICO case
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 11.0 KB
Line 
1! -*- Mode: f90 -*-
2PROGRAM interpol
3!!!
4!!! Generates interpolation weights from source (src) grid to destination (dst) grid
5!!! Most of the configuration is done by the xml file
6!!!
7   USE xios
8   USE mod_wait
9   USE mpi
10   IMPLICIT NONE
11   ! INCLUDE "mpif.h"
12   INTEGER, PARAMETER :: rl = SELECTED_REAL_KIND(15,307) !< Default real precision (double)
13   INTEGER :: ierr, comm, rank
14
15   TYPE (xios_context) :: ctx_hdl
16   INTEGER :: ni_src, nj_src !< Dimensions of the source grid
17   INTEGER :: ni_dst, nj_dst !< Dimensions of the source grid
18   REAL (kind=rl), ALLOCATABLE :: imask_src (:,:), imask_dst (:,:), area_src(:,:), area_dst(:,:)
19   REAL (kind=rl), ALLOCATABLE :: lon_src(:,:), lat_src(:,:), field_src(:,:)
20   LOGICAL, ALLOCATABLE :: lmask_src (:,:), lmask_dst (:,:)
21   INTEGER :: nout = 0, jf
22   CHARACTER (LEN=20) :: nchar, type_src, type_dst
23   LOGICAL :: l_mask_src = .TRUE. , l_mask_dst = .TRUE., l_use_area = .FALSE.
24   LOGICAL :: l_src_is2D = .FALSE., l_src_is1D = .FALSE., l_src_is0D = .FALSE., l_dst_is2D = .FALSE. , l_dst_is1D = .FALSE., l_dst_is0D = .FALSE.
25   !
26   REAL (kind=rl), PARAMETER :: rpi = ACOS ( -1.0_rl)
27   REAL (kind=rl), PARAMETER :: rad = rpi / 180.0_rl
28   !
29   INTEGER :: narg, ja !< To handle command line arguments
30   CHARACTER (LEN=80) :: cmd_arg
31   !
32   !! SVN information
33   CHARACTER (LEN = LEN ( "$Author$" )) :: &
34      &   SVN_Author   =  "$Author$"
35   CHARACTER (LEN = LEN ( "$Date$" )) :: &
36      &   SVN_Date     =  "$Date$"
37   CHARACTER (LEN = LEN ( "$Revision$")) :: &
38      &   SVN_Revision =  "$Revision$"
39   CHARACTER (LEN = LEN ( "$Id$" )) :: &
40      &   SVN_Id       =  "$Id$"
41   CHARACTER (LEN = LEN ( & ! Strange line arrangement to match the max width (132 chars)
42    & "$HeadURL$" &
43    & )) :: &
44    & SVN_HeadURL  =  &
45    & "$HeadURL$"
46
47   !!
48   !< Initialisation
49   WRITE (UNIT=nout, FMT="('-- start interpol')" ) 
50   CALL MPI_INIT (ierr)
51
52   CALL init_wait
53   CALL xios_initialize ("interpol", return_comm=comm)
54
55   CALL MPI_COMM_RANK (comm, rank, ierr)
56   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi_comm_rank')" ) rank
57   !CALL MPI_COMM_SIZE (comm, size, ierr)
58
59   !< Read command line
60   narg = iargc ()
61   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- narg : ', 1I3)" ) rank, narg
62
63   !< Read arguments on command line
64   DO ja = 1, narg
65      CALL getarg ( ja, cmd_arg )
66      SELECT CASE ( TRIM (cmd_arg) )
67      CASE ( '--mask_src=true'  ) ; l_mask_src = .TRUE.
68      CASE ( '--mask_src=false' ) ; l_mask_src = .FALSE.
69      CASE ( '--mask_dst=true'  ) ; l_mask_dst = .TRUE.
70      CASE ( '--mask_dst=false' ) ; l_mask_dst = .FALSE.
71      CASE ( '--use_area=true'  ) ; l_use_area = .TRUE.
72      CASE ( '--use_area=no'    ) ; l_use_area = .FALSE.
73      END SELECT
74   END DO
75
76   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mask_src : ', 1L)" ) rank, l_mask_src
77   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mask_dst : ', 1L)" ) rank, l_mask_dst
78   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- use_area : ', 1L)" ) rank, l_use_area
79
80   !< Context interpol_read : read masks
81   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Starting interpol_read')" ) rank
82   CALL xios_context_initialize  ("interpol_read", comm)
83   CALL xios_get_handle          ("interpol_read", ctx_hdl)
84   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Set current context interpol_read')" ) rank
85   CALL xios_set_current_context (ctx_hdl)
86   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Close context definition interpol_read')" ) rank
87   CALL xios_close_context_definition ()
88
89   !< Read characteristics of the source grid
90   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading domain_src')" ) rank
91   CALL xios_get_domain_attr ("domain_src", ni=ni_src, nj=nj_src)
92   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading domain_src ', 6I9)") rank, ni_src, nj_src
93   ALLOCATE ( lon_src (ni_src, nj_src), lat_src (ni_src, nj_src) )
94   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading lon lat src ', 6I9)") rank, SIZE (lon_src), &
95      &        SHAPE (lon_src), SIZE (lat_src), SHAPE (lat_src)
96   !!
97   !
98   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading type src ', 6I9)") rank
99   CALL xios_get_domain_attr ("domain_src", TYPE=type_src )
100   SELECT CASE ( TRIM (type_src))
101   CASE ( "rectilinear" )
102      l_src_is2D = .TRUE.
103      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading lon lat src rectilinear', 6I9)") rank
104      CALL xios_get_domain_attr ("domain_src", lonvalue_1d=lon_src(:,1), latvalue_1d=lat_src(1,:) )
105      lon_src (:,:) = SPREAD ( lon_src(:,1), DIM=2, ncopies=nj_src)
106      lat_src (:,:) = SPREAD ( lat_src(1,:), DIM=1, ncopies=nj_src)
107   CASE default
108      IF ( nj_src == 1 .AND. ni_src == 1) THEN
109         l_src_is0D = .TRUE.
110         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading lon lat src 0D ', 6I9)") rank
111         CALL xios_get_domain_attr ("domain_src", lonvalue_2d=lon_src(:,:), latvalue_2d=lat_src(:,:) )
112      ELSE IF ( nj_src == 1 ) THEN
113         l_src_is1D = .TRUE.
114         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading lon lat src 1D ', 6I9)") rank
115         CALL xios_get_domain_attr ("domain_src", lonvalue_1d=lon_src(:,1), latvalue_1d=lat_src(:,1) )
116      ELSE
117         l_src_is2D = .TRUE.
118         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading lon lat src 2D ', 6I9)") rank
119         CALL xios_get_domain_attr ("domain_src", lonvalue_2d=lon_src(:,:), latvalue_2d=lat_src(:,:) )
120      ENDIF
121   END SELECT
122   
123   !< Read mask on the source grid
124   ALLOCATE ( imask_src (ni_src, nj_src), lmask_src (ni_src, nj_src) )
125   IF ( l_mask_src ) THEN
126      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Receive field mask_src')" ) rank
127      CALL xios_recv_field ("mask_src", imask_src)
128      lmask_src = .FALSE.
129      WHERE (imask_src > 0.5 ) lmask_src = .TRUE.
130   ELSE
131      imask_src (:,:) = 1 ; lmask_src (:,:) = .TRUE.
132   ENDIF
133   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Counting mask_src : ', 1I8)" ) rank, COUNT(lmask_src)
134
135   !< Read area on the source grid
136   IF ( l_use_area ) THEN
137      IF ( l_src_is0D ) ALLOCATE ( area_src (  1   , 1     ) )
138      IF ( l_src_is1D ) ALLOCATE ( area_src (  1   , ni_src) )
139      IF ( l_src_is2D ) ALLOCATE ( area_src (ni_src, nj_src) )
140      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Receive field area_src')" ) rank
141      CALL xios_recv_field ("area_src", area_src(:,:))
142      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Sum area_src : ', 2E15.3)" ) rank, SUM(area_src)
143   ENDIF
144   
145   !< Read characteristics of the destination grid
146   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading domain_dst')" ) rank
147   CALL xios_get_domain_attr ("domain_dst", ni=ni_dst, nj=nj_dst)
148   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading domain_dst ', 6I7)") rank, ni_dst, nj_dst
149
150   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Reading type dst ', 6I9)") rank
151   CALL xios_get_domain_attr ("domain_dst", TYPE=type_dst )
152   SELECT CASE ( TRIM (type_dst))
153   CASE ( "rectilinear" )
154      l_dst_is2D = .TRUE.
155      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Rectilinear dest domain 2D')") rank
156   CASE default
157      IF ( nj_dst == 1 .AND. ni_dst == 1) THEN
158         l_dst_is0D = .TRUE.
159         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Domain dest 0D')") rank
160      ELSE IF ( nj_dst == 1 ) THEN
161         l_dst_is1D = .TRUE.
162         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Domain dest 1D')") rank
163      ELSE
164         l_dst_is2D = .TRUE.
165         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Domain dest 2D')") rank
166      ENDIF
167   END SELECT
168   
169   !< Read mask on the destination grid
170   ALLOCATE ( imask_dst (ni_dst, nj_dst), lmask_dst (ni_dst, nj_dst) )
171   IF ( l_mask_dst ) THEN
172      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Receive field mask_dst')" ) rank 
173      CALL xios_recv_field ("mask_dst", imask_dst)
174      lmask_dst = .FALSE.
175      WHERE (imask_dst > 0.5) lmask_dst = .TRUE.
176   ELSE
177      imask_dst (:,:) = 1 ; lmask_dst (:,:) = .TRUE.
178   ENDIF
179   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- counting mask_dst : ', 1I8)" ) rank, COUNT(lmask_dst)
180
181   !< Read area on the destination grid
182   IF ( l_use_area ) THEN
183      IF ( l_dst_is0D ) ALLOCATE ( area_dst (  1   , 1     ) )
184      IF ( l_dst_is1D ) ALLOCATE ( area_dst (  1   , ni_dst) )
185      IF ( l_dst_is2D ) ALLOCATE ( area_dst (ni_dst, nj_dst) )
186      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- receive field area_dst')" ) rank 
187      CALL xios_recv_field ("area_dst", area_dst(:,:))
188      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Sum area_dst : ', 2E15.3)" ) rank, SUM(area_dst)
189   ENDIF
190   
191   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- End interpol_read (1)')") rank
192   CALL xios_context_finalize ()
193   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- End interpol_read (2)')") rank
194   
195   !< Context interpol run : generates weights, interpolate mask from source to destination
196   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Starting interpol_run')" ) rank
197   CALL xios_context_initialize  ("interpol_run", comm)
198   CALL xios_get_handle          ("interpol_run", ctx_hdl)
199   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Set context interpol_run')" ) rank
200   CALL xios_set_current_context (ctx_hdl)
201
202   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Get attributes domain_src')" ) rank
203   IF ( l_mask_src ) CALL xios_set_domain_attr ("domain_src", mask_2d=lmask_src)
204   IF ( l_use_area ) CALL xios_set_domain_attr ("domain_src", area   =area_src )
205   !CALL xios_set_domain_attr ("domain_src", radius= 6371229.0 )
206   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Get attributes domain_dst')" ) rank
207   IF ( l_mask_dst ) CALL xios_set_domain_attr ("domain_dst", mask_2d=lmask_dst)
208   IF ( l_use_area ) CALL xios_set_domain_attr ("domain_dst", area   =area_dst )
209   !CALL xios_set_domain_attr ("domain_dst", radius= 6371229.0 )
210   
211   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Close context definition')" ) rank
212   CALL xios_close_context_definition ()
213
214   CALL xios_update_calendar (1)
215   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Send field mask_src')" ) rank
216   CALL xios_send_field ("mask_src", imask_src)
217
218   !!< Creates analytic fields and interpolate
219   ALLOCATE ( field_src (ni_src, nj_src) )
220
221   DO jf = 1, 6
222      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- Working on test case ', 1I2.2)" ) rank, jf
223      SELECT CASE ( jf)
224      CASE ( 1) ; field_src (:,:) = REAL ( imask_src, kind=rl)
225      CASE ( 2) ; field_src (:,:) = SIN ( rad * lat_src(:,:) )
226      CASE ( 3) ; field_src (:,:) = COS ( rad * lat_src(:,:) )
227      CASE ( 4) ; field_src (:,:) = COS ( rad * lat_src(:,:) ) * COS ( rad * lon_src(:,:) )
228      CASE ( 5) ; field_src (:,:) = SIN ( rad * lon_src(:,:) )
229      CASE ( 6) ; field_src (:,:) = SIN ( rad * lon_src(:,:) + rpi/2.0_rl )
230      END SELECT
231      WRITE (UNIT=nchar, FMT='("field", 1I2.2, "_src")' ) jf
232      CALL xios_send_field ( TRIM(nchar), field_src )
233   END DO
234   
235   !!<
236   
237   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- End interpol_run')" ) rank
238   CALL xios_context_finalize ()
239
240   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi_comm_free')" ) rank
241   CALL MPI_COMM_FREE (comm, ierr)
242
243   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- xios finalize')" ) rank
244   CALL xios_finalize ()
245
246   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi finalize')" ) rank
247   CALL MPI_FINALIZE (ierr)
248
249   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- The end')" ) rank
250
251END PROGRAM interpol
252
253
254
255
256
Note: See TracBrowser for help on using the repository browser.