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

Last change on this file since 4153 was 4153, checked in by omamce, 3 years ago

O.M. : include areas correction

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