source: trunk/WRF.COMMON/WRFV2/external/io_esmf/ext_esmf_read_field.F90 @ 2756

Last change on this file since 2756 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 17.6 KB
Line 
1
2!$$$here...  TBH:  remove duplication between ext_esmf_read_field and
3!$$$here...  TBH:  ext_esmf_write_field
4
5!$$$here...  TBH:  how to deal with time?  (via current ESMF_Clock)
6!$$$here...  TBH:  to begin, use it as an error check! 
7
8
9!--- read_field
10SUBROUTINE ext_esmf_read_field ( DataHandle , DateStr , VarName , Field , FieldType , Comm , IOComm, &
11                                 DomainDesc , MemoryOrder , Stagger , DimNames ,              &
12                                 DomainStart , DomainEnd ,                                    &
13                                 MemoryStart , MemoryEnd ,                                    &
14                                 PatchStart , PatchEnd ,                                      &
15                                 Status )
16  USE module_ext_esmf
17  IMPLICIT NONE
18  INTEGER       ,INTENT(IN)    :: DataHandle
19  CHARACTER*(*) ,intent(inout) :: DateStr
20  CHARACTER*(*) ,intent(inout) :: VarName
21  integer       ,intent(inout) :: FieldType
22  integer       ,intent(inout) :: Comm
23  integer       ,intent(inout) :: IOComm
24  integer       ,intent(inout) :: DomainDesc
25  character*(*) ,intent(inout) :: MemoryOrder
26  character*(*) ,intent(inout) :: Stagger
27  character*(*) ,intent(inout) :: DimNames(*)
28  integer       ,intent(inout) :: DomainStart(*), DomainEnd(*)
29  integer       ,intent(inout) :: MemoryStart(*), MemoryEnd(*)
30  integer       ,intent(inout) :: PatchStart(*),  PatchEnd(*)
31  REAL          ,INTENT(INOUT) :: Field(*)
32  integer       ,intent(out)   :: Status
33  ! Local declarations
34  INTEGER :: ids,ide,jds,jde,kds,kde
35  INTEGER :: ims,ime,jms,jme,kms,kme
36  INTEGER :: ips,ipe,jps,jpe,kps,kpe
37  TYPE(ESMF_State), POINTER :: importstate
38  TYPE(ESMF_Field) :: tmpField
39  TYPE(ESMF_Array) :: tmpArray
40  TYPE(ESMF_ArraySpec) :: arrayspec
41  TYPE(ESMF_DataKind) :: esmf_kind
42  TYPE(ESMF_DataType) :: esmf_type
43  TYPE(ESMF_RelLoc) :: horzRelloc
44  REAL(ESMF_KIND_R4), POINTER :: data_esmf_real_ptr(:,:)
45  REAL(ESMF_KIND_R4), POINTER :: tmp_esmf_r4_ptr(:,:)
46  INTEGER(ESMF_KIND_I4), POINTER :: data_esmf_int_ptr(:,:)
47  INTEGER, PARAMETER :: esmf_rank = 2
48  INTEGER :: DomainEndFull(esmf_rank), idefull, jdefull, ict, i, j
49  INTEGER :: PatchEndFull(esmf_rank), ipefull, jpefull
50  ! esmf_counts is redundant.  remove it as soon as ESMF_ArrayCreate no
51  ! longer requires it
52  INTEGER :: esmf_counts(esmf_rank)
53  INTEGER :: rc
54  LOGICAL, EXTERNAL :: has_char
55  character*256 mess
56!$$$DEBUG
57INTEGER, SAVE :: numtimes=0   ! track number of calls
58CHARACTER(LEN=256) :: timestamp
59!REAL :: debug_real(MemoryStart(1):MemoryEnd(1),MemoryStart(2):MemoryEnd(2))
60!$$$END DEBUG
61
62  IF ( .NOT. int_valid_handle( DataHandle ) ) THEN
63    CALL wrf_error_fatal("ext_esmf_read_field: invalid data handle" )
64  ENDIF
65  IF ( .NOT. int_handle_in_use( DataHandle ) ) THEN
66    CALL wrf_error_fatal("ext_esmf_read_field: DataHandle not opened" )
67  ENDIF
68  IF ( .NOT. opened_for_read( DataHandle ) ) THEN
69    CALL wrf_error_fatal("ext_esmf_read_field: DataHandle not opened for read" )
70  ENDIF
71
72write(mess,*)'ext_esmf_read_field ',DataHandle, TRIM(DateStr), TRIM(VarName)
73call wrf_debug( 300, TRIM(mess) )
74
75  IF      ( FieldType .EQ. WRF_REAL ) THEN
76    esmf_type = ESMF_DATA_REAL
77    esmf_kind = ESMF_R4
78  ELSE IF ( FieldType .EQ. WRF_DOUBLE ) THEN
79!    esmf_type = ESMF_DATA_REAL
80!    esmf_kind = ESMF_R8
81    CALL wrf_error_fatal( 'ext_esmf_read_field, WRF_DOUBLE not yet supported')
82  ELSE IF ( FieldType .EQ. WRF_INTEGER ) THEN
83    esmf_type = ESMF_DATA_INTEGER
84    esmf_kind = ESMF_I4
85!$$$ implement this (below)
86    CALL wrf_error_fatal( 'ext_esmf_read_field, WRF_INTEGER not yet implemented')
87  ELSE IF ( FieldType .EQ. WRF_LOGICAL ) THEN
88    CALL wrf_error_fatal( 'ext_esmf_read_field, WRF_LOGICAL not yet supported')
89  ENDIF
90
91  ims = MemoryStart(1) ; ime = MemoryEnd(1)
92  jms = MemoryStart(2) ; jme = MemoryEnd(2)
93  kms = MemoryStart(3) ; kme = MemoryEnd(3)
94
95  ips = PatchStart(1) ; ipe = PatchEnd(1)
96  jps = PatchStart(2) ; jpe = PatchEnd(2)
97  kps = PatchStart(3) ; kpe = PatchEnd(3)
98
99  ids = DomainStart(1) ; ide = DomainEnd(1)
100  jds = DomainStart(2) ; jde = DomainEnd(2)
101  kds = DomainStart(3) ; kde = DomainEnd(3)
102
103  ! For now, treat all arrays as 2D... 
104!$$$ Eventually, use ../io_netcdf subroutines Transpose() and reorder()
105!$$$ (and etc.) to handle general array ranks and index orderings. 
106!$$$ Some copies of these exist in ../../frame/module_io.F. 
107!$$$ Then use ESMF_ArrayDataMap class to handle index mapping. 
108  IF ( kms /= kme ) THEN
109    CALL wrf_error_fatal( 'ext_esmf_read_field:  rank > 2 not yet supported')
110  ENDIF
111
112! The non-staggered variables come in at one-less than
113! domain dimensions, but io_esmf is currently hacked to use full
114! domain spec, so adjust if not staggered.
115! $$$ TBD:  Remove EndFull hackery once ESMF can support staggered
116! $$$ TBD:  grids in regional models.  (This hack works around the current
117! $$$ TBD:  need to use only larger staggered dimensions for ESMF_Arrays.) 
118  CALL ioesmf_endfullhack( esmf_rank, DomainEnd, PatchEnd, Stagger, &
119                           DomainEndFull, PatchEndFull )
120  idefull = DomainEndFull(1)
121  jdefull = DomainEndFull(2)
122  ipefull = PatchEndFull(1)
123  jpefull = PatchEndFull(2)
124
125write(mess,*) ' ext_esmf_read_field: okay_to_read: ', DataHandle, okay_to_read(DataHandle)
126call wrf_debug( 300, TRIM(mess) )
127
128  ! case 1: the file is opened for read but not committed ("training")
129  IF ( .NOT. okay_to_read( DataHandle ) )  THEN
130
131    ! Training:  build the ESMF import state
132write(mess,*) ' ext_esmf_read_field: TRAINING READ:  DataHandle = ', DataHandle
133call wrf_debug( 300, TRIM(mess) )
134
135    ! First, build the ESMF_Grid for this DataHandle, if it does not
136    ! already exist
137    CALL ioesmf_create_grid( DataHandle, esmf_rank, MemoryOrder, Stagger,      &
138                             DomainStart(1:esmf_rank), DomainEnd(1:esmf_rank), &
139                             MemoryStart(1:esmf_rank), MemoryEnd(1:esmf_rank), &
140                             PatchStart(1:esmf_rank), PatchEnd(1:esmf_rank) )
141    ! Grab the current importState and add to it...
142    CALL ESMF_ImportStateGetCurrent( importstate, rc )
143    IF ( rc /= ESMF_SUCCESS ) THEN
144      CALL wrf_error_fatal("ext_esmf_read_field, training:  ESMF_ImportStateGetCurrent failed" )
145    ENDIF
146! BEGIN DOESNOTWORK
147! The following code does not work for reasons as-yet unknown. 
148! A likely suspect is lbounds and ubounds which fail in other interfaces in
149! ESMF 2.2.0rp1 ... 
150    ! Build ESMF objects... 
151    ! Build an ESMF_ArraySpec.  The use of ESMF_ArraySpec and ESMF_Array
152    ! objects allows some of the code that follows to be type-kind-independent. 
153!    CALL ESMF_ArraySpecSet(arrayspec, rank=esmf_rank, type=esmf_type, &
154!                                      kind=esmf_kind, rc=rc)
155!    IF ( rc /= ESMF_SUCCESS ) THEN
156!      CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_ArraySpecSet failed" )
157!    ENDIF
158    ! Build an ESMF_Array
159    ! Implementation note:  since we do not yet have full control over how
160    ! ESMF chooses to lay out a "patch" within "memory", we must copy by
161    ! hand.  (Reasons include lack of support in ESMF for asymmetric halos,
162    ! addition of "extra" rows/columns to optimize alignment on some machines,
163    ! handling of periodic boundary conditions, etc.)  Thus, there
164    ! is no point in using larger "memory" sizes to build the array -- patch
165    ! is fine.  Also, since we must copy anyway, might as well let ESMF manage
166    ! the memory for simplicity. 
167!$$$ Once ESMF can match WRF memory-patch mapping, replace this with a more
168!$$$ efficient solution that does not require a copy. 
169! $$$ esmf_counts is redundant.  Remove it as soon as ESMF_ArrayCreate no
170! $$$ longer requires it. 
171!    esmf_counts(1:esmf_rank) = DomainEndFull(1:esmf_rank) - &
172!                               DomainStart(1:esmf_rank) + 1
173!    tmpArray = ESMF_ArrayCreate(arrayspec, counts=esmf_counts,      &
174!                                lbounds=DomainStart(1:esmf_rank),   &
175!                                ubounds=DomainEndFull(1:esmf_rank), &
176!                                rc=rc)
177!    IF ( rc /= ESMF_SUCCESS ) THEN
178!      WRITE(mess,*) ' ext_esmf_read_field: ESMF_ArrayCreate failed, rc = ', rc
179!      CALL wrf_error_fatal( TRIM(mess) )
180!    ENDIF
181    ! Determine grid staggering for this Field
182!    IF ( has_char( Stagger, 'x' ) .AND. has_char( Stagger, 'y' ) ) THEN
183!      CALL wrf_error_fatal( &
184!        "ext_esmf_read_field:  ESMF does not yet support XY staggering for C-grid" )
185!    ELSE IF ( has_char( Stagger, 'x' ) ) THEN
186!      horzrelloc=ESMF_CELL_WFACE
187!    ELSE IF ( has_char( Stagger, 'y' ) ) THEN
188!      horzrelloc=ESMF_CELL_SFACE
189!    ELSE
190!      horzrelloc=ESMF_CELL_CENTER
191!    ENDIF
192    ! Build an ESMF_Field
193    ! Note:  though it is counter-intuitive, ESMF uses
194    ! shallow-copy-masquerading-as-reference to implement the
195    ! pseudo-equivalent of POINTER assignment under-the-hood.  What this means
196    ! here is that it is OK to pass deep object tmpArray into
197    ! ESMF_FieldCreate() and then return from this subroutine.  Even though
198    ! tmpArray goes out of scope, it is OK.  However, if tmpArray were to be
199    ! modified after this call, the changes would not be guaranteed to always
200    ! appear in tmpField.  It works that way now, but ESMF Core team has
201    ! plans that may make it break in the future.  Build-it, attach-it,
202    ! flush-it will work.  Build-it, attach-it, modify-it, flush-it may not
203    ! always work. 
204    ! "Pie, pie and a fox..." 
205    ! Note:  unique Field name is required by ESMF_StateAddField(). 
206!$$$here...  use CF "standard_name" once the WRF Registry supports it
207!    tmpField = ESMF_FieldCreate( grid( DataHandle )%ptr, tmpArray,          &
208!                                 copyflag=ESMF_DATA_REF,                    &
209!                                 horzrelloc=horzrelloc, name=TRIM(VarName), &
210!                                 rc=rc )
211! END DOESNOTWORK
212    !$$$here...  This is a complete HACK for debugging!!  Need to compute
213    !$$$here...  horzrelloc from Stagger as above... 
214    horzrelloc=ESMF_CELL_CENTER
215    !$$$ TODO:  Add code for other data types here... 
216    ALLOCATE( tmp_esmf_r4_ptr(ips:ipefull,jps:jpefull) )
217    CALL wrf_debug ( 100, 'ext_esmf_read_field: calling ESMF_FieldCreate' )
218    tmpField = ESMF_FieldCreate(         &
219                 grid( DataHandle )%ptr, &
220                 tmp_esmf_r4_ptr,        &
221                 copyflag=ESMF_DATA_REF, &
222                 horzrelloc=horzrelloc,  &
223                 name=TRIM(VarName),     &
224                 rc=rc )
225    IF ( rc /= ESMF_SUCCESS ) THEN
226      WRITE(mess,*) ' ext_esmf_read_field: ESMF_FieldCreate failed, rc = ', rc
227      CALL wrf_error_fatal( TRIM(mess) )
228    ENDIF
229    CALL wrf_debug ( 100, 'ext_esmf_read_field: back from ESMF_FieldCreate' )
230    WRITE(mess,*) 'ext_esmf_read_field: tmp_esmf_r4_ptr(',         &
231      LBOUND(tmp_esmf_r4_ptr,1),':',UBOUND(tmp_esmf_r4_ptr,1),',', &
232      LBOUND(tmp_esmf_r4_ptr,2),':',UBOUND(tmp_esmf_r4_ptr,2),')'
233    CALL wrf_debug ( 100 , TRIM(mess) )
234    ! Add the Field to the import state... 
235!$$$here...  for now, just build ESMF_Fields and stuff them in
236!$$$here...  later, use a single ESMF_Bundle
237    CALL ESMF_StateAddField( importstate, tmpField, rc=rc )
238    IF ( rc /= ESMF_SUCCESS ) THEN
239      CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_StateAddField failed" )
240    ENDIF
241write(mess,*) ' ext_esmf_read_field: END TRAINING READ:  DataHandle = ', DataHandle
242call wrf_debug( 300, TRIM(mess) )
243
244  ! case 2: opened for read and committed
245  ELSE IF ( okay_to_read( DataHandle ) )  THEN
246
247write(mess,*) ' ext_esmf_read_field: ACTUAL READ:  DataHandle = ', DataHandle
248call wrf_debug( 300, TRIM(mess) )
249!$$$DEBUG
250! count calls...
251numtimes = numtimes + 1
252CALL get_current_time_string( timestamp )
253!$$$END DEBUG
254
255    ! read:  extract data from the ESMF import state
256    ! Grab the current importState
257    CALL ESMF_ImportStateGetCurrent( importstate, rc )
258    IF ( rc /= ESMF_SUCCESS ) THEN
259      CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_ImportStateGetCurrent failed" )
260    ENDIF
261    ! grab the Field
262    CALL ESMF_StateGetField( importstate, fieldName=TRIM(VarName), &
263                             field=tmpfield, rc=rc )
264    IF ( rc /= ESMF_SUCCESS ) THEN
265      CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_StateGetField failed" )
266    ENDIF
267!$$$DEBUG
268CALL wrf_debug ( 100, 'ext_esmf_read_field '//TRIM(VarName)//':  calling ESMF_FieldPrint( tmpField ) 1' )
269CALL ESMF_FieldPrint( tmpField, rc=rc )
270CALL wrf_debug ( 100, 'ext_esmf_read_field '//TRIM(VarName)//':  back from ESMF_FieldPrint( tmpField ) 1' )
271!$$$END DEBUG
272
273    ! grab a pointer to the import state data and copy data into Field
274    IF      ( FieldType .EQ. WRF_REAL ) THEN
275      CALL ESMF_FieldGetDataPointer( tmpField, data_esmf_real_ptr, &
276                                     ESMF_DATA_REF, rc=rc )
277      IF ( rc /= ESMF_SUCCESS ) THEN
278        CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_FieldGetDataPointer(r4) failed" )
279      ENDIF
280      IF ( ( PatchStart(1)   /= LBOUND(data_esmf_real_ptr,1) ) .OR. &
281           ( PatchEndFull(1) /= UBOUND(data_esmf_real_ptr,1) ) .OR. &
282           ( PatchStart(2)   /= LBOUND(data_esmf_real_ptr,2) ) .OR. &
283           ( PatchEndFull(2) /= UBOUND(data_esmf_real_ptr,2) ) ) THEN
284        WRITE( mess,* ) 'ESMF_FieldGetDataPointer bounds mismatch',          &
285          __FILE__ ,                                                         &
286          ', line ',                                                         &
287          __LINE__ ,                                                         &
288          ', ips:ipe,jps:jpe = ',PatchStart(1),':',PatchEndFull(1),',',      &
289                                 PatchStart(2),':',PatchEndFull(2),          &
290          ', data_esmf_real_ptr(BOUNDS) = ',                                 &
291          LBOUND(data_esmf_real_ptr,1),':',UBOUND(data_esmf_real_ptr,1),',', &
292          LBOUND(data_esmf_real_ptr,2),':',UBOUND(data_esmf_real_ptr,2)
293        CALL wrf_error_fatal ( TRIM(mess) )
294      ENDIF
295!$$$DEBUG
296WRITE( mess,* ) 'DEBUG:  ext_esmf_read_field:  ips:ipe,jps:jpe = ',  &
297  ips,':',ipe,',',jps,':',jpe,                                       &
298  ', data_esmf_real_ptr(BOUNDS) = ',                                 &
299  LBOUND(data_esmf_real_ptr,1),':',UBOUND(data_esmf_real_ptr,1),',', &
300  LBOUND(data_esmf_real_ptr,2),':',UBOUND(data_esmf_real_ptr,2)
301CALL wrf_debug( 300, TRIM(mess) )
302!DO j= jms, jme
303!  DO i= ims, ime
304!    debug_real(i,j) = -(i*1000.0 + j)/100000.0     ! obvious bad value for debugging
305!  ENDDO
306!ENDDO
307!debug_real(ips:ipe,jps:jpe) = data_esmf_real_ptr(ips:ipe,jps:jpe)
308!CALL wrf_debug( 100, 'DEBUG:  ext_esmf_read_field:  writing DEBUG1_WRFcmp_import'//TRIM(VarName)//'_'//TRIM(timestamp) )
309!OPEN( UNIT=985, FILE='DEBUG1_WRFcmp_import'//TRIM(VarName)//'_'//TRIM(timestamp), FORM='formatted' )
310!WRITE (985,'(a,a,i4)') TRIM(VarName),' ',numtimes
311!DO j = jps, jpe
312!  DO i = ips, ipe
313!    WRITE (985,*) '(',i,',',j,'):  ',debug_real(i,j)
314!  ENDDO
315!ENDDO
316!CLOSE (985)
317!$$$END DEBUG
318      CALL ioesmf_extract_data_real( data_esmf_real_ptr, Field,            &
319                                     ips, ipefull, jps, jpefull, kps, kpe, &
320                                     ims, ime, jms, jme, kms, kme )
321!$$$DEBUG
322!ict = 0
323!DO j= jms, jme
324!  DO i= ims, ime
325!    ict = ict + 1
326!    IF ( (i<ips) .OR. (i>ipe) .OR. (j<jps) .OR. (j>jpe) ) THEN
327!      debug_real(i,j) = -(i*1000.0 + j)/100000.0     ! obvious bad value for debugging
328!    ELSE
329!      debug_real(i,j) = Field(ict)
330!    ENDIF
331!  ENDDO
332!ENDDO
333!CALL wrf_debug( 100, 'DEBUG:  ext_esmf_read_field:  writing DEBUG1_WRFcmp_read_Field'//TRIM(VarName)//'_'//TRIM(timestamp) )
334!OPEN( UNIT=985, FILE='DEBUG1_WRFcmp_read_Field'//TRIM(VarName)//'_'//TRIM(timestamp), FORM='formatted' )
335!WRITE (985,'(a,a,i4)') TRIM(VarName),' ',numtimes
336!DO j = jps, jpe
337!  DO i = ips, ipe
338!    WRITE (985,*) '(',i,',',j,'):  ',debug_real(i,j)
339!  ENDDO
340!ENDDO
341!CLOSE (985)
342!$$$END DEBUG
343    ELSE IF ( FieldType .EQ. WRF_INTEGER ) THEN
344      CALL ESMF_FieldGetDataPointer( tmpField, data_esmf_int_ptr, &
345                                     ESMF_DATA_REF, rc=rc )
346      IF ( rc /= ESMF_SUCCESS ) THEN
347        CALL wrf_error_fatal("ext_esmf_read_field:  ESMF_FieldGetDataPointer(i4) failed" )
348      ENDIF
349      IF ( ( PatchStart(1)   /= LBOUND(data_esmf_int_ptr,1) ) .OR. &
350           ( PatchEndFull(1) /= UBOUND(data_esmf_int_ptr,1) ) .OR. &
351           ( PatchStart(2)   /= LBOUND(data_esmf_int_ptr,2) ) .OR. &
352           ( PatchEndFull(2) /= UBOUND(data_esmf_int_ptr,2) ) ) THEN
353        WRITE( mess,* ) 'ESMF_FieldGetDataPointer bounds mismatch',        &
354          __FILE__ ,                                                       &
355          ', line ',                                                       &
356          __LINE__ ,                                                       &
357          ', ips:ipe,jps:jpe = ',PatchStart(1),':',PatchEndFull(1),',',    &
358                                 PatchStart(2),':',PatchEndFull(2),        &
359          ', data_esmf_int_ptr(BOUNDS) = ',                                &
360          LBOUND(data_esmf_int_ptr,1),':',UBOUND(data_esmf_int_ptr,1),',', &
361          LBOUND(data_esmf_int_ptr,2),':',UBOUND(data_esmf_int_ptr,2)
362        CALL wrf_error_fatal ( TRIM(mess) )
363      ENDIF
364      CALL ioesmf_extract_data_int( data_esmf_int_ptr, Field,             &
365                                    ips, ipefull, jps, jpefull, kps, kpe, &
366                                    ims, ime, jms, jme, kms, kme )
367    ENDIF
368write(mess,*) ' ext_esmf_read_field: END ACTUAL READ:  DataHandle = ', DataHandle
369call wrf_debug( 300, TRIM(mess) )
370
371  ENDIF
372
373!$$$DEBUG
374CALL wrf_debug ( 100, 'ext_esmf_read_field '//TRIM(VarName)//':  calling ESMF_FieldPrint( tmpField ) 2' )
375CALL ESMF_FieldPrint( tmpField, rc=rc )
376CALL wrf_debug ( 100, 'ext_esmf_read_field '//TRIM(VarName)//':  back from ESMF_FieldPrint( tmpField ) 2' )
377!$$$END DEBUG
378
379  Status = 0
380
381  RETURN
382
383END SUBROUTINE ext_esmf_read_field
384
Note: See TracBrowser for help on using the repository browser.