2 !$$$here... TBH: remove duplication between ext_esmf_read_field and
3 !$$$here... TBH: ext_esmf_write_field
5 !$$$here... TBH: how to deal with time? (via current ESMF_Clock)
6 !$$$here... TBH: to begin, use it as an error check!
10 SUBROUTINE ext_esmf_write_field ( DataHandle , DateStr , VarName , Field , FieldType , Comm , IOComm, &
11 DomainDesc , MemoryOrder , Stagger , DimNames , &
12 DomainStart , DomainEnd , &
13 MemoryStart , MemoryEnd , &
14 PatchStart , PatchEnd , &
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
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 :: exportstate
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
52 INTEGER :: esmf_counts(esmf_rank)
54 LOGICAL, EXTERNAL :: has_char
57 INTEGER, SAVE :: numtimes=0 ! track number of calls
58 CHARACTER(LEN=256) :: timestamp
59 !REAL :: debug_real(MemoryStart(1):MemoryEnd(1),MemoryStart(2):MemoryEnd(2))
62 IF ( .NOT. int_valid_handle( DataHandle ) ) THEN
63 CALL wrf_error_fatal("ext_esmf_write_field: invalid data handle" )
65 IF ( .NOT. int_handle_in_use( DataHandle ) ) THEN
66 CALL wrf_error_fatal("ext_esmf_write_field: DataHandle not opened" )
68 IF ( .NOT. opened_for_write( DataHandle ) ) THEN
69 CALL wrf_error_fatal("ext_esmf_write_field: DataHandle not opened for write" )
72 write(mess,*)'ext_esmf_write_field ',DataHandle, TRIM(DateStr), TRIM(VarName)
73 call wrf_debug( 300, TRIM(mess) )
75 IF ( FieldType .EQ. WRF_REAL ) THEN
76 esmf_type = ESMF_DATA_REAL
78 ELSE IF ( FieldType .EQ. WRF_DOUBLE ) THEN
79 ! esmf_type = ESMF_DATA_REAL
81 CALL wrf_error_fatal( 'ext_esmf_write_field, WRF_DOUBLE not yet supported')
82 ELSE IF ( FieldType .EQ. WRF_INTEGER ) THEN
83 esmf_type = ESMF_DATA_INTEGER
85 !$$$ implement this (below)
86 CALL wrf_error_fatal( 'ext_esmf_write_field, WRF_INTEGER not yet implemented')
87 ELSE IF ( FieldType .EQ. WRF_LOGICAL ) THEN
88 CALL wrf_error_fatal( 'ext_esmf_write_field, WRF_LOGICAL not yet supported')
91 ims = MemoryStart(1) ; ime = MemoryEnd(1)
92 jms = MemoryStart(2) ; jme = MemoryEnd(2)
93 kms = MemoryStart(3) ; kme = MemoryEnd(3)
95 ips = PatchStart(1) ; ipe = PatchEnd(1)
96 jps = PatchStart(2) ; jpe = PatchEnd(2)
97 kps = PatchStart(3) ; kpe = PatchEnd(3)
99 ids = DomainStart(1) ; ide = DomainEnd(1)
100 jds = DomainStart(2) ; jde = DomainEnd(2)
101 kds = DomainStart(3) ; kde = DomainEnd(3)
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_write_field: rank > 2 not yet supported')
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)
125 write(mess,*) ' ext_esmf_write_field: okay_to_write: ', DataHandle, okay_to_write(DataHandle)
126 call wrf_debug( 300, TRIM(mess) )
128 ! case 1: the file is opened for write but not committed ("training")
129 IF ( .NOT. okay_to_write( DataHandle ) ) THEN
131 ! Training: build the ESMF export state
132 write(mess,*) ' ext_esmf_write_field: TRAINING WRITE: DataHandle = ', DataHandle
133 call wrf_debug( 300, TRIM(mess) )
135 ! First, build the ESMF_Grid for this DataHandle, if it does not
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 exportState and add to it...
142 CALL ESMF_ExportStateGetCurrent( exportstate, rc )
143 IF ( rc /= ESMF_SUCCESS ) THEN
144 CALL wrf_error_fatal("ext_esmf_write_field, training: ESMF_ExportStateGetCurrent failed" )
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
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_write_field: ESMF_ArraySpecSet failed" )
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), &
177 ! IF ( rc /= ESMF_SUCCESS ) THEN
178 ! WRITE(mess,*) ' ext_esmf_write_field: ESMF_ArrayCreate failed, rc = ', rc
179 ! CALL wrf_error_fatal( TRIM(mess) )
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_write_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
190 ! horzrelloc=ESMF_CELL_CENTER
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
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), &
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_write_field: calling ESMF_FieldCreate' )
218 tmpField = ESMF_FieldCreate( &
219 grid( DataHandle )%ptr, &
221 copyflag=ESMF_DATA_REF, &
222 horzrelloc=horzrelloc, &
223 name=TRIM(VarName), &
225 IF ( rc /= ESMF_SUCCESS ) THEN
226 WRITE(mess,*) ' ext_esmf_write_field: ESMF_FieldCreate failed, rc = ', rc
227 CALL wrf_error_fatal( TRIM(mess) )
229 CALL wrf_debug ( 100, 'ext_esmf_write_field: back from ESMF_FieldCreate' )
230 WRITE(mess,*) 'ext_esmf_write_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 export 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( exportstate, tmpField, rc=rc )
238 IF ( rc /= ESMF_SUCCESS ) THEN
239 CALL wrf_error_fatal("ext_esmf_write_field: ESMF_StateAddField failed" )
241 write(mess,*) ' ext_esmf_write_field: END TRAINING WRITE: DataHandle = ', DataHandle
242 call wrf_debug( 300, TRIM(mess) )
244 ! case 2: opened for write and committed
245 ELSE IF ( okay_to_write( DataHandle ) ) THEN
247 write(mess,*) ' ext_esmf_write_field: ACTUAL WRITE: DataHandle = ', DataHandle
248 call wrf_debug( 300, TRIM(mess) )
251 numtimes = numtimes + 1
252 CALL get_current_time_string( timestamp )
255 ! write: insert data into the ESMF export state
256 ! Grab the current exportState
257 CALL ESMF_ExportStateGetCurrent( exportstate, rc )
258 IF ( rc /= ESMF_SUCCESS ) THEN
259 CALL wrf_error_fatal("ext_esmf_write_field: ESMF_ExportStateGetCurrent failed" )
262 CALL ESMF_StateGetField( exportstate, fieldName=TRIM(VarName), &
263 field=tmpfield, rc=rc )
264 IF ( rc /= ESMF_SUCCESS ) THEN
265 CALL wrf_error_fatal("ext_esmf_write_field: ESMF_StateGetField failed" )
268 CALL wrf_debug ( 100, 'ext_esmf_write_field '//TRIM(VarName)//': calling ESMF_FieldPrint( tmpField ) 1' )
269 CALL ESMF_FieldPrint( tmpField, rc=rc )
270 CALL wrf_debug ( 100, 'ext_esmf_write_field '//TRIM(VarName)//': back from ESMF_FieldPrint( tmpField ) 1' )
273 ! grab a pointer to the export state data and copy data from 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_write_field: ESMF_FieldGetDataPointer(r4) failed" )
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', &
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) )
296 WRITE( mess,* ) 'DEBUG: ext_esmf_write_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)
301 CALL wrf_debug( 300, TRIM(mess) )
306 ! IF ( (i<ips) .OR. (i>ipe) .OR. (j<jps) .OR. (j>jpe) ) THEN
307 ! debug_real(i,j) = -(i*1000.0 + j)/100000.0 ! obvious bad value for debugging
309 ! debug_real(i,j) = Field(ict)
313 !CALL wrf_debug( 100, 'DEBUG: ext_esmf_write_field: writing DEBUG1_WRFcmp_write_Field'//TRIM(VarName)//'_'//TRIM(timestamp) )
314 !OPEN( UNIT=985, FILE='DEBUG1_WRFcmp_write_Field'//TRIM(VarName)//'_'//TRIM(timestamp), FORM='formatted' )
315 !WRITE (985,'(a,a,i4)') TRIM(VarName),' ',numtimes
318 ! WRITE (985,*) '(',i,',',j,'): ',debug_real(i,j)
323 CALL ioesmf_insert_data_real( Field, data_esmf_real_ptr, &
324 ips, ipefull, jps, jpefull, kps, kpe, &
325 ims, ime, jms, jme, kms, kme )
329 ! debug_real(i,j) = -(i*1000.0 + j)/100000.0 ! obvious bad value for debugging
332 !debug_real(ips:ipe,jps:jpe) = data_esmf_real_ptr(ips:ipe,jps:jpe)
333 !CALL wrf_debug( 100, 'DEBUG: ext_esmf_write_field: writing DEBUG1_WRFcmp_export'//TRIM(VarName)//'_'//TRIM(timestamp) )
334 !OPEN( UNIT=985, FILE='DEBUG1_WRFcmp_export'//TRIM(VarName)//'_'//TRIM(timestamp), FORM='formatted' )
335 !WRITE (985,'(a,a,i4)') TRIM(VarName),' ',numtimes
338 ! WRITE (985,*) '(',i,',',j,'): ',debug_real(i,j)
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_write_field: ESMF_FieldGetDataPointer(i4) failed" )
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', &
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) )
364 CALL ioesmf_insert_data_int( Field, data_esmf_int_ptr, &
365 ips, ipefull, jps, jpefull, kps, kpe, &
366 ims, ime, jms, jme, kms, kme )
368 write(mess,*) ' ext_esmf_write_field: END ACTUAL WRITE: DataHandle = ', DataHandle
369 call wrf_debug( 300, TRIM(mess) )
374 CALL wrf_debug ( 100, 'ext_esmf_write_field '//TRIM(VarName)//': calling ESMF_FieldPrint( tmpField ) 2' )
375 CALL ESMF_FieldPrint( tmpField, rc=rc )
376 CALL wrf_debug ( 100, 'ext_esmf_write_field '//TRIM(VarName)//': back from ESMF_FieldPrint( tmpField ) 2' )
383 END SUBROUTINE ext_esmf_write_field