Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hco_geotools_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCO\_GeoTools\_Mod contains a collection of
9 : ! helper routines for extracting geographical information. These
10 : ! routines are based upon GEOS-5 data and may need to be revised
11 : ! for other met. fields!
12 : !\\
13 : !\\
14 : ! !INTERFACE:
15 : !
16 : MODULE HCO_GeoTools_Mod
17 : !
18 : ! !USES:
19 : !
20 : USE HCO_Error_Mod
21 :
22 : IMPLICIT NONE
23 : PRIVATE
24 : !
25 : ! !PUBLIC MEMBER FUNCTIONS:
26 : !
27 : PUBLIC :: HCO_LandType
28 : PUBLIC :: HCO_ValidateLon
29 : PUBLIC :: HCO_GetSUNCOS
30 : PUBLIC :: HCO_GetHorzIJIndex
31 : PUBLIC :: HCO_CalcVertGrid
32 : PUBLIC :: HCO_SetPBLm
33 : ! PUBLIC :: HCO_CalcPBLlev
34 :
35 : INTERFACE HCO_LandType
36 : MODULE PROCEDURE HCO_LandType_Dp
37 : MODULE PROCEDURE HCO_LandType_Sp
38 : END INTERFACE HCO_LandType
39 :
40 : INTERFACE HCO_ValidateLon
41 : MODULE PROCEDURE HCO_ValidateLon_Dp
42 : MODULE PROCEDURE HCO_ValidateLon_Sp
43 : END INTERFACE HCO_ValidateLon
44 :
45 : ! INTERFACE HCO_CalcPBLlev
46 : ! MODULE PROCEDURE HCO_CalcPBLlev2D
47 : ! MODULE PROCEDURE HCO_CalcPBLlev3D
48 : ! END INTERFACE HCO_CalcPBLlev
49 : !
50 : ! !PRIVATE MEMBER FUNCTIONS:
51 : !
52 : PRIVATE:: HCO_LandType_Dp
53 : PRIVATE:: HCO_LandType_Sp
54 : PRIVATE:: HCO_ValidateLon_Dp
55 : PRIVATE:: HCO_ValidateLon_Sp
56 : !
57 : ! !REVISION HISTORY:
58 : ! 18 Dec 2013 - C. Keller - Initialization
59 : ! See https://github.com/geoschem/hemco for complete history
60 : !EOP
61 : !------------------------------------------------------------------------------
62 : !BOC
63 : CONTAINS
64 : !EOC
65 : !------------------------------------------------------------------------------
66 : ! Harmonized Emissions Component (HEMCO) !
67 : !------------------------------------------------------------------------------
68 : !BOP
69 : !
70 : ! !IROUTINE: HCO_LandType_Sp
71 : !
72 : ! !DESCRIPTION: Function HCO\_LANDTYPE returns the land type based upon
73 : ! GMAO land type fractions. Result is 0=water, 1=land, 2=ice, where ice includes
74 : ! over both ocean and land, and land includes lakes. Inputs are in single precision.
75 : !\\
76 : !\\
77 : ! !INTERFACE:
78 : !
79 0 : FUNCTION HCO_LandType_Sp( FRLAND, FRLANDIC, FROCEAN, FRSEAICE, FRLAKE ) &
80 : Result ( LandType )
81 : !
82 : ! !INPUT PARAMETERS:
83 : !
84 : REAL(sp), INTENT(IN) :: FRLAND ! Fraction of grid-box covered in land
85 : REAL(sp), INTENT(IN) :: FRLANDIC ! Fraction of grid-box covered in land ice
86 : REAL(sp), INTENT(IN) :: FROCEAN ! Fraction of grid-box covered in ocean
87 : REAL(sp), INTENT(IN) :: FRSEAICE ! Fraction of grid-box covered in sea ice
88 : REAL(sp), INTENT(IN) :: FRLAKE ! Fraction of grid-box covered in lake
89 : !
90 : ! !RETURN VALUE
91 : !
92 : INTEGER :: LandType ! Land type: 0=ocean, 1=land, 2=ice (ocn,land)
93 : !
94 : ! !REMARKS:
95 : ! Land type index is based on GMAO field LWI, with modification
96 : ! to classify land ice as ice.
97 : !
98 : ! !REVISION HISTORY:
99 : ! 18 Dec 2013 - C. Keller - Initialization!
100 : ! See https://github.com/geoschem/hemco for complete history
101 : !EOP
102 : !------------------------------------------------------------------------------
103 : !BOC
104 : !
105 : ! !DEFINED PARAMETERS:
106 : !
107 : ! Threshold at which a grid-box is considered ice
108 : REAL(sp), PARAMETER :: frac_classify_land_ice_as_ice = 0.5_sp
109 :
110 : !--------------------------
111 : ! HCO_LANDTYPE begins here
112 : !--------------------------
113 :
114 : ! Start with the surface type category definitions based on the GMAO
115 : ! definition for land-water-ice index (LWI), which is 0=ocean (non-ice),
116 : ! 1=land (includes lakes and ice), 2=ice (over ocean only). Qualify
117 : ! land type as location of maximum fraction.
118 : LandType = MAXLOC( (/ FRLAND + FRLANDIC + FRLAKE, &
119 : FRSEAICE, &
120 0 : FROCEAN - FRSEAICE /), 1 )
121 0 : IF ( LandType == 3 ) LandType = 0
122 :
123 : ! Change land type to ice if sufficient ice over land
124 0 : IF ( FRLANDIC >= frac_classify_land_ice_as_ice ) THEN
125 0 : LandType = 2
126 : ENDIF
127 :
128 0 : END FUNCTION HCO_LandType_Sp
129 : !EOC
130 : !------------------------------------------------------------------------------
131 : ! Harmonized Emissions Component (HEMCO) !
132 : !------------------------------------------------------------------------------
133 : !BOP
134 : !
135 : ! !IROUTINE: HCO_LandType_Dp
136 : !
137 : ! !DESCRIPTION: Function HCO\_LANDTYPE returns the land type based upon
138 : ! GMAO land type fractions. Result is 0=water, 1=land, 2=ice, where ice includes
139 : ! over both ocean and land, and land includes lakes. Inputs are in double precision.
140 : !\\
141 : !\\
142 : ! !INTERFACE:
143 : !
144 0 : FUNCTION HCO_LandType_Dp( FRLAND, FRLANDIC, FROCEAN, FRSEAICE, FRLAKE ) &
145 : Result ( LandType )
146 : !
147 : ! !INPUT PARAMETERS:
148 : !
149 : REAL(dp), INTENT(IN) :: FRLANDIC ! Fraction of grid-box covered in land ice
150 : REAL(dp), INTENT(IN) :: FRLAND ! Fraction of grid-box covered in land
151 : REAL(dp), INTENT(IN) :: FROCEAN ! Fraction of grid-box covered in ocean
152 : REAL(dp), INTENT(IN) :: FRSEAICE ! Fraction of grid-box covered in sea ice
153 : REAL(dp), INTENT(IN) :: FRLAKE ! Fraction of grid-box covered in lake
154 : !
155 : ! !RETURN VALUE:
156 : !
157 : INTEGER :: LandType ! Land type: 0=ocean, 1=land, 2=ice (ocn,land)
158 : !
159 : ! !REMARKS:
160 : ! Land type index is based on GMAO field LWI, with modification
161 : ! to classify land ice as ice.
162 : !
163 : ! !REVISION HISTORY:
164 : ! 18 Dec 2013 - C. Keller - Initialization
165 : ! See https://github.com/geoschem/hemco for complete history
166 : !EOP
167 : !------------------------------------------------------------------------------
168 : !BOC
169 : !
170 : ! !DEFINED PARAMETERS::
171 : !
172 : ! Threshold at which a grid-box with land ice is considered ice
173 : REAL(dp), PARAMETER :: frac_classify_land_ice_as_ice = 0.5_dp
174 :
175 : !--------------------------
176 : ! HCO_LANDTYPE begins here
177 : !--------------------------
178 :
179 : ! Start with the surface type category definitions based on the GMAO
180 : ! definition for land-water-ice index (LWI), which is 0=ocean (non-ice),
181 : ! 1=land (includes lakes and ice), 2=ice (over ocean only). Qualify
182 : ! land type as location of maximum fraction.
183 : LandType = MAXLOC( (/ FRLAND + FRLANDIC + FRLAKE, &
184 : FRSEAICE, &
185 0 : FROCEAN - FRSEAICE /), 1 )
186 0 : IF ( LandType == 3 ) LandType = 0
187 :
188 : ! Change land type to ice if sufficient ice over land
189 0 : IF ( FRLANDIC >= frac_classify_land_ice_as_ice ) THEN
190 0 : LandType = 2
191 : ENDIF
192 :
193 0 : END FUNCTION HCO_LandType_Dp
194 : !EOC
195 : !------------------------------------------------------------------------------
196 : ! Harmonized Emissions Component (HEMCO) !
197 : !------------------------------------------------------------------------------
198 : !BOP
199 : !
200 : ! !IROUTINE: HCO_ValidateLon_Sp
201 : !
202 : ! !DESCRIPTION: Subroutine HCO\_ValidateLon\_Sp ensures that the passed
203 : ! single precision longitude axis LON is steadily increasing.
204 : !\\
205 : !\\
206 : ! !INTERFACE:
207 : !
208 0 : SUBROUTINE HCO_ValidateLon_Sp ( HcoState, NLON, LON, RC )
209 : !
210 : ! !USES:
211 : !
212 : USE HCO_STATE_MOD, ONLY : HCO_STATE
213 : !
214 : ! !INPUT/OUTPUT PARAMETERS:
215 : !
216 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
217 : INTEGER, INTENT(IN ) :: NLON ! # of lons
218 : !
219 : ! !INPUT/OUTPUT PARAMETERS:
220 : !
221 : REAL(sp), INTENT(INOUT) :: LON(NLON) ! longitude axis
222 : INTEGER, INTENT(INOUT) :: RC ! Return code
223 : !
224 : ! !REVISION HISTORY:
225 : ! 16 Jul 2014 - C. Keller - Initialization
226 : ! See https://github.com/geoschem/hemco for complete history
227 : !EOP
228 : !------------------------------------------------------------------------------
229 : !BOC
230 : !
231 : ! !LOCAL VARIABLES:
232 : !
233 : INTEGER :: I, CNT
234 : LOGICAL :: REDO
235 : INTEGER, PARAMETER :: MAXIT = 10
236 :
237 : !--------------------------
238 : ! HCO_ValidateLon_Sp begins here
239 : !--------------------------
240 :
241 0 : REDO = .TRUE.
242 0 : CNT = 0
243 :
244 0 : DO WHILE ( REDO )
245 :
246 : ! Exit w/ error after 10 iterations
247 0 : CNT = CNT + 1
248 0 : IF ( CNT > MAXIT ) THEN
249 : CALL HCO_ERROR ( '>10 iterations', RC, &
250 0 : THISLOC='HCO_ValidateLon (HCO_GEOTOOLS_MOD.F90)' )
251 0 : RETURN
252 : ENDIF
253 :
254 0 : DO I = 1, NLON
255 :
256 : ! If we reach the last grid box, all values are steadily
257 : ! increasing (otherwise, the loop would have been exited).
258 0 : IF ( I == NLON ) THEN
259 : REDO = .FALSE.
260 : EXIT
261 : ENDIF
262 :
263 : ! Check if next lon value is lower, in which case we subtract
264 : ! a value of 360 (degrees) from all longitude values up to
265 : ! this point. Then repeat the lookup (from the beginning).
266 0 : IF ( LON(I+1) < LON(I) ) THEN
267 0 : LON(1:I) = LON(1:I) - 360.0_sp
268 : EXIT
269 : ENDIF
270 :
271 : ENDDO !I
272 : ENDDO ! REDO
273 :
274 : ! Leave w/ success
275 0 : RC = HCO_SUCCESS
276 :
277 : END SUBROUTINE HCO_ValidateLon_Sp
278 : !EOC
279 : !------------------------------------------------------------------------------
280 : ! Harmonized Emissions Component (HEMCO) !
281 : !------------------------------------------------------------------------------
282 : !BOP
283 : !
284 : ! !IROUTINE: HCO_ValidateLon_Dp
285 : !
286 : ! !DESCRIPTION: Subroutine HCO\_ValidateLon\_Sp ensures that the passed
287 : ! double precision longitude axis LON is steadily increasing.
288 : !\\
289 : !\\
290 : ! !INTERFACE:
291 : !
292 0 : SUBROUTINE HCO_ValidateLon_Dp ( HcoState, NLON, LON, RC )
293 : !
294 : ! !USES:
295 : !
296 : USE HCO_STATE_MOD, ONLY : HCO_STATE
297 : !
298 : ! !INPUT/OUTPUT PARAMETERS:
299 : !
300 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
301 : INTEGER, INTENT(IN ) :: NLON ! # of lons
302 : !
303 : !
304 : ! !INPUT/OUTPUT PARAMETERS:
305 : !
306 : REAL(dp), INTENT(INOUT) :: LON(NLON) ! longitude axis
307 : INTEGER, INTENT(INOUT) :: RC ! Return code
308 : !
309 : ! !REVISION HISTORY:
310 : ! 16 Jul 2014 - C. Keller - Initialization
311 : ! See https://github.com/geoschem/hemco for complete history
312 : !EOP
313 : !------------------------------------------------------------------------------
314 : !BOC
315 : INTEGER :: I, CNT
316 : LOGICAL :: REDO
317 : INTEGER, PARAMETER :: MAXIT = 10
318 :
319 : !--------------------------
320 : ! HCO_ValidateLon_Dp begins here
321 : !--------------------------
322 :
323 0 : REDO = .TRUE.
324 0 : CNT = 0
325 :
326 0 : DO WHILE ( REDO )
327 :
328 : ! Exit w/ error after 10 iterations
329 0 : CNT = CNT + 1
330 0 : IF ( CNT > MAXIT ) THEN
331 : CALL HCO_ERROR ( '>10 iterations', RC, &
332 0 : THISLOC='HCO_ValidateLon (HCO_GEOTOOLS_MOD.F90)' )
333 0 : RETURN
334 : ENDIF
335 :
336 0 : DO I = 1, NLON
337 :
338 : ! If we reach the last grid box, all values are steadily
339 : ! increasing (otherwise, the loop would have been exited).
340 0 : IF ( I == NLON ) THEN
341 : REDO = .FALSE.
342 : EXIT
343 : ENDIF
344 :
345 : ! Check if next lon value is lower, in which case we subtract
346 : ! a value of 360 (degrees) from all longitude values up to
347 : ! this point. Then repeat the lookup (from the beginning).
348 0 : IF ( LON(I+1) < LON(I) ) THEN
349 0 : LON(1:I) = LON(1:I) - 360.0_dp
350 : EXIT
351 : ENDIF
352 :
353 : ENDDO !I
354 : ENDDO ! REDO
355 :
356 : ! Leave w/ success
357 0 : RC = HCO_SUCCESS
358 :
359 : END SUBROUTINE HCO_ValidateLon_Dp
360 : !EOC
361 : !------------------------------------------------------------------------------
362 : ! Harmonized Emissions Component (HEMCO) !
363 : !------------------------------------------------------------------------------
364 : !BOP
365 : !
366 : ! !IROUTINE: HCO_GetSUNCOS
367 : !
368 : ! !DESCRIPTION: Subroutine HCO\_GetSUNCOS calculates the solar zenith angle
369 : ! for the given date.
370 : !\\
371 : !\\
372 : ! !INTERFACE:
373 : !
374 0 : SUBROUTINE HCO_GetSUNCOS( HcoState, SUNCOS, DT, RC )
375 : !
376 : ! !USES
377 : !
378 : USE HCO_STATE_MOD, ONLY : HCO_STATE
379 : USE HCO_CLOCK_MOD, ONLY : HcoClock_Get
380 : USE HCO_CLOCK_MOD, ONLY : HcoClock_GetLocal
381 : !
382 : ! !INPUT/OUTPUT PARAMETERS:
383 : !
384 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
385 : INTEGER, INTENT(IN ) :: DT ! Time shift relative to current date [hrs]
386 : !
387 : ! !OUTPUT PARAMETERS:
388 : !
389 : REAL(hp), INTENT( OUT) :: SUNCOS(HcoState%NX,HcoState%NY)
390 : !
391 : ! !INPUT/OUTPUT PARAMETERS:
392 : !
393 : INTEGER, INTENT(INOUT) :: RC ! Return code
394 : !
395 : ! !REVISION HISTORY:
396 : ! 22 May 2015 - C. Keller - Initial version, based on GEOS-Chem's dao_mod.F.
397 : ! See https://github.com/geoschem/hemco for complete history
398 : !EOP
399 : !------------------------------------------------------------------------------
400 : !BOC
401 : !
402 : ! !LOCAL VARIABLES:
403 : !
404 : INTEGER :: I, J, DOY, HOUR
405 : LOGICAL :: ERR
406 : REAL(hp) :: YMID_R, S_YMID_R, C_YMID_R
407 : REAL(hp) :: R, DEC
408 : REAL(hp) :: S_DEC, C_DEC
409 : REAL(hp) :: SC, LHR
410 : REAL(hp) :: AHR
411 :
412 : ! Coefficients for solar declination angle
413 : REAL(hp), PARAMETER :: A0 = 0.006918e+0_hp
414 : REAL(hp), PARAMETER :: A1 = 0.399912e+0_hp
415 : REAL(hp), PARAMETER :: A2 = 0.006758e+0_hp
416 : REAL(hp), PARAMETER :: A3 = 0.002697e+0_hp
417 : REAL(hp), PARAMETER :: B1 = 0.070257e+0_hp
418 : REAL(hp), PARAMETER :: B2 = 0.000907e+0_hp
419 : REAL(hp), PARAMETER :: B3 = 0.000148e+0_hp
420 :
421 : CHARACTER(LEN=255) :: LOC
422 :
423 : !-------------------------------
424 : ! HCO_GetSUNCOS starts here!
425 : !-------------------------------
426 0 : LOC = 'HCO_GetSUNCOS (HCO_GEOTOOLS_MOD.F90)'
427 :
428 : ! Get current time information
429 0 : CALL HcoClock_Get( HcoState%Clock, cDOY=DOY, cH=HOUR, RC=RC )
430 0 : IF ( RC /= HCO_SUCCESS ) THEN
431 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
432 0 : RETURN
433 : ENDIF
434 :
435 : ! Add time adjustment
436 0 : HOUR = HOUR + DT
437 :
438 : ! Make sure HOUR is within valid range (0-24)
439 0 : IF ( HOUR < 0 ) THEN
440 0 : HOUR = HOUR + 24
441 0 : DOY = DOY - 1
442 0 : ELSEIF ( HOUR > 23 ) THEN
443 0 : HOUR = HOUR - 24
444 0 : DOY = DOY + 1
445 : ENDIF
446 :
447 : ! Make sure DOY is within valid range of 1 to 365
448 0 : DOY = MAX(MIN(DOY,365),1)
449 :
450 : ! Path length of earth's orbit traversed since Jan 1 [radians]
451 0 : R = ( 2e+0_hp * HcoState%Phys%PI / 365e+0_hp ) * DBLE( DOY - 1 )
452 :
453 : ! Solar declination angle (low precision formula) [radians]
454 : DEC = A0 - A1*COS( R ) + B1*SIN( R ) &
455 : - A2*COS( 2e+0_hp*R ) + B2*SIN( 2e+0_hp*R ) &
456 0 : - A3*COS( 3e+0_hp*R ) + B3*SIN( 3e+0_hp*R )
457 :
458 : ! Pre-compute sin & cos of DEC outside of DO loops (for efficiency)
459 0 : S_DEC = SIN( DEC )
460 0 : C_DEC = COS( DEC )
461 :
462 : ! Init
463 0 : ERR = .FALSE.
464 :
465 : !=================================================================
466 : ! Compute cosine of solar zenith angle
467 : !=================================================================
468 : !$OMP PARALLEL DO &
469 : !$OMP DEFAULT( SHARED ) &
470 : !$OMP PRIVATE( I, J, YMID_R, S_YMID_R, C_YMID_R ) &
471 : !$OMP PRIVATE( LHR, AHR, SC, RC )
472 0 : DO J = 1, HcoState%NY
473 0 : DO I = 1, HcoState%NX
474 :
475 : ! Latitude of grid box [radians]
476 0 : YMID_R = HcoState%Grid%YMID%Val(I,J) * HcoState%Phys%PI_180
477 :
478 : ! Pre-compute sin & cos of DEC outside of I loop (for efficiency)
479 0 : S_YMID_R = SIN( YMID_R )
480 0 : C_YMID_R = COS( YMID_R )
481 :
482 : !==============================================================
483 : ! Compute cosine of SZA at the midpoint of the chem timestep
484 : ! Required for photolysis, chemistry, emissions, drydep
485 : !==============================================================
486 :
487 : !-----------------------------------------------------------------------------
488 : ! Prior to 3/2/17:
489 : ! Seb Eastham suggested to comment out the call to HcoClock_GetLocal. If
490 : ! the Voronoi timezones are used, this will compute the timezones on political
491 : ! boundaries and not strictly on longitude. This will cause funny results.
492 : ! Replace this with a strict longitudinal local time. (bmy, 3/27/17)
493 : ! ! Local time [hours] at box (I,J) at the midpt of the chem timestep
494 : ! CALL HcoClock_GetLocal ( HcoState, I, J, cH=LHR, RC=RC )
495 : !
496 : ! IF ( RC /= HCO_SUCCESS ) THEN
497 : ! ERR = .TRUE.
498 : ! EXIT
499 : ! ENDIF
500 : !-----------------------------------------------------------------------------
501 : ! Prior to 3/2/17:
502 : ! Seb Eastham says that HOUR (in the new formula below) already contains DT.
503 : ! so we need to comment this out and just use HOUR + LONGITUDE/15.
504 : ! (bmy, 3/2/17)
505 : ! ! Adjust for time shift
506 : ! LHR = LHR + DT
507 : !----------------------------------------------------------------------------
508 :
509 : ! Compute local time as UTC + longitude/15 (bmy, 3/2/17)
510 0 : LHR = HOUR + ( HcoState%Grid%XMid%Val(I,J) / 15.0_hp )
511 :
512 0 : IF ( LHR < 0.0_hp ) LHR = LHR + 24.0_hp
513 0 : IF ( LHR >= 24.0_hp ) LHR = LHR - 24.0_hp
514 :
515 : ! Hour angle at box (I,J) [radians]
516 0 : AHR = ABS( LHR - 12.0_hp ) * 15.0_hp * HcoState%Phys%PI_180
517 :
518 : ! Corresponding cosine( SZA ) at box (I,J) [unitless]
519 : SC = ( S_YMID_R * S_DEC ) &
520 0 : + ( C_YMID_R * C_DEC * COS( AHR ) )
521 :
522 : ! COS(SZA) at the current time
523 0 : SUNCOS(I,J) = SC
524 :
525 : ENDDO
526 : ENDDO
527 : !$OMP END PARALLEL DO
528 :
529 : ! Check error status
530 : IF ( ERR ) THEN
531 : CALL HCO_ERROR ( &
532 : 'Cannot calculate SZA', RC, &
533 : THISLOC='HCO_GetSUNCOS (hco_geotools_mod.F90)' )
534 : RETURN
535 : ENDIF
536 :
537 : ! Leave w/ success
538 0 : RC = HCO_SUCCESS
539 :
540 : END SUBROUTINE HCO_GetSUNCOS
541 : !EOC
542 : #if defined(ESMF_)
543 : !------------------------------------------------------------------------------
544 : ! Harmonized Emissions Component (HEMCO) !
545 : !------------------------------------------------------------------------------
546 : !BOP
547 : !
548 : ! !IROUTINE: HCO_GetHorzIJIndex
549 : !
550 : ! !DESCRIPTION: Function HCO\_GetHorzIJIndex returns the grid box index for
551 : ! the given longitude (deg E, -180...180), and latitude (deg N, -90...90).
552 : !\\
553 : !\\
554 : ! !INTERFACE:
555 : !
556 : SUBROUTINE HCO_GetHorzIJIndex( HcoState, N, Lon, Lat, idx, jdx, RC )
557 : !
558 : ! !USES
559 : !
560 : #include "MAPL_Generic.h"
561 : USE ESMF
562 : USE MAPLBase_Mod
563 : USE HCO_STATE_MOD, ONLY : HCO_STATE
564 : !
565 : ! !INPUT PARAMETERS:
566 : !
567 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
568 : INTEGER, INTENT(IN ) :: N
569 : REAL(hp), INTENT(IN ) :: Lon(N)
570 : REAL(hp), INTENT(IN ) :: Lat(N)
571 : !
572 : ! !INPUT/OUTPUT PARAMETERS:
573 : !
574 : INTEGER, INTENT(INOUT) :: RC
575 : !
576 : ! !OUTPUT PARAMETERS:
577 : !
578 : INTEGER, INTENT( OUT) :: IDX(N), JDX(N)
579 : !
580 : ! !REVISION HISTORY:
581 : ! 04 Jun 2015 - C. Keller - Initial version
582 : ! See https://github.com/geoschem/hemco for complete history
583 : !EOP
584 : !------------------------------------------------------------------------------
585 : !BOC
586 : !
587 : ! !LOCAL VARIABLES:
588 : !
589 : REAL :: LonR(N), LatR(N)
590 : REAL, PARAMETER :: radToDeg = 57.2957795
591 : TYPE(ESMF_Grid) :: Grid
592 :
593 : ! Defined Iam and STATUS
594 : __Iam__("HCO_GetHorzIJIndex (hco_geotools_mod.F90)")
595 :
596 : !-------------------------------
597 : ! HCO_GetHorzIJIndex begins here
598 : !-------------------------------
599 :
600 : ! Get grid
601 : ASSERT_(ASSOCIATED(HcoState%GridComp))
602 : CALL ESMF_GridCompGet( HcoState%GridComp, Grid=Grid, __RC__ )
603 :
604 : ! Shadow variables
605 : LonR(:) = Lon / radToDeg
606 : LatR(:) = Lat / radToDeg
607 :
608 : ! Get indeces
609 : CALL MAPL_GetHorzIJIndex( npts=N, II=idx, JJ=jdx, &
610 : lon=LonR, lat=LatR, Grid=Grid, &
611 : __RC__)
612 :
613 : !!! old version of MAPL:
614 : ! CALL MAPL_GetHorzIJIndex(N,idx,jdx,LonR,LatR,Grid=Grid,__RC__)
615 :
616 : ! Return w/ success
617 : RC = HCO_SUCCESS
618 :
619 : END SUBROUTINE HCO_GetHorzIJIndex
620 : !EOC
621 : #else
622 : !------------------------------------------------------------------------------
623 : ! Harmonized Emissions Component (HEMCO) !
624 : !------------------------------------------------------------------------------
625 : !BOP
626 : !
627 : ! !IROUTINE: HCO_GetHorzIJIndex
628 : !
629 : ! !DESCRIPTION: Function HCO\_GetHorzIJIndex returns the grid box index for
630 : ! the given longitude (deg E, -180...180), and latitude (deg N, -90...90).
631 : !\\
632 : !\\
633 : ! !INTERFACE:
634 : !
635 0 : SUBROUTINE HCO_GetHorzIJIndex( HcoState, N, Lon, Lat, idx, jdx, RC )
636 : !
637 : ! !USES:
638 : !
639 : USE HCO_STATE_MOD, ONLY : HCO_STATE
640 : !
641 : ! !INPUT PARAMETERS:
642 : !
643 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
644 : INTEGER, INTENT(IN ) :: N
645 : REAL(hp), INTENT(IN ) :: Lon(N)
646 : REAL(hp), INTENT(IN ) :: Lat(N)
647 : !
648 : ! !INPUT/OUTPUT PARAMETERS:
649 : !
650 : INTEGER, INTENT(INOUT) :: RC
651 : !
652 : ! !OUTPUT PARAMETERS:
653 : !
654 : INTEGER, INTENT( OUT) :: IDX(N), JDX(N)
655 : !
656 : ! !REVISION HISTORY:
657 : ! 04 Jun 2015 - C. Keller - Initial version
658 : ! See https://github.com/geoschem/hemco for complete history
659 : !EOP
660 : !------------------------------------------------------------------------------
661 : !BOC
662 : !
663 : ! !LOCAL VARIABLES:
664 : !
665 : INTEGER :: I, J, L, FOUND
666 : REAL(hp) :: iLon1, iLon2
667 : REAL(hp) :: iLat1, iLat2
668 : REAL(hp) :: delta
669 :
670 : !-------------------------------
671 : ! HCO_GetHorzIJIndex begins here
672 : !-------------------------------
673 :
674 : ! Initialize
675 0 : IDX(:) = -1
676 0 : JDX(:) = -1
677 0 : FOUND = 0
678 :
679 : ! do for every grid box
680 0 : DO J = 1, HcoState%NY
681 0 : DO I = 1, HcoState%NX
682 :
683 : ! Get grid edges for this box
684 :
685 : ! Longitude edges
686 0 : IF ( ASSOCIATED(HcoState%Grid%XEDGE%Val) ) THEN
687 0 : iLon1 = HcoState%Grid%XEDGE%Val(I, J)
688 0 : iLon2 = HcoState%Grid%XEDGE%Val(I+1,J)
689 :
690 : ! Approximate from mid points
691 : ELSE
692 0 : iLon1 = HcoState%Grid%XMID%Val(I,J)
693 0 : IF ( I < HcoState%NX ) THEN
694 0 : delta = HcoState%Grid%XMID%Val(I+1,J)-iLon1
695 : ELSE
696 0 : delta = iLon1 - HcoState%Grid%XMID%Val(I-1,J)
697 : ENDIF
698 0 : iLon2 = iLon1 + (delta / 2.0_hp)
699 0 : iLon1 = iLon1 - (delta / 2.0_hp)
700 : ENDIF
701 :
702 : ! Latitude edges
703 0 : IF ( ASSOCIATED(HcoState%Grid%YEDGE%Val) ) THEN
704 0 : iLat1 = HcoState%Grid%YEDGE%Val(I,J)
705 0 : iLat2 = HcoState%Grid%YEDGE%Val(I,J+1)
706 :
707 : ! Approximate from mid points
708 : ELSE
709 0 : iLat1 = HcoState%Grid%YMID%Val(I,J)
710 0 : IF ( J < HcoState%NY ) THEN
711 0 : delta = HcoState%Grid%YMID%Val(I,J+1)-iLat1
712 : ELSE
713 0 : delta = iLat1 - HcoState%Grid%YMID%Val(I,J-1)
714 : ENDIF
715 0 : iLat2 = iLat1 + (delta / 2.0_hp)
716 0 : iLat1 = iLat1 - (delta / 2.0_hp)
717 : ENDIF
718 :
719 : ! Check if it's within this box
720 0 : DO L = 1, N
721 0 : IF ( IDX(L) > 0 ) CYCLE
722 :
723 0 : IF ( Lon(L) >= HcoState%Grid%XEDGE%Val(I, J ) .AND. &
724 0 : Lon(L) <= HcoState%Grid%XEDGE%Val(I+1,J ) .AND. &
725 0 : Lat(L) >= HcoState%Grid%YEDGE%Val(I ,J ) .AND. &
726 0 : Lat(L) <= HcoState%Grid%YEDGE%Val(I ,J+1) ) THEN
727 0 : IDX(L) = I
728 0 : JDX(L) = J
729 0 : FOUND = FOUND + 1
730 0 : IF ( FOUND == N ) EXIT
731 : ENDIF
732 : ENDDO
733 : ENDDO
734 : ENDDO
735 :
736 : ! Return w/ success
737 0 : RC = HCO_SUCCESS
738 :
739 0 : END SUBROUTINE HCO_GetHorzIJIndex
740 : !EOC
741 : #endif
742 : !------------------------------------------------------------------------------
743 : ! Harmonized Emissions Component (HEMCO) !
744 : !------------------------------------------------------------------------------
745 : !BOP
746 : !
747 : ! !IROUTINE: HCO_CalcVertGrid
748 : !
749 : ! !DESCRIPTION: Function HCO\_CalcVertGrid calculates the vertical grid
750 : ! quantities surface pressure PSFC [Pa], surface geopotential height ZSFC
751 : ! [m], grid box height BXHEIGHT [m], and pressure edges PEDGE [Pa]. Any of
752 : ! these fields can be passed explicitly to the routine, in which case these
753 : ! fields are being used. If not passed through the routine (i.e. if the
754 : ! corresponding input argument pointer is nullified), the field is searched
755 : ! in the HEMCO configuration file. If not found in the configuration file,
756 : ! the field is approximated from other quantities (if possible). For example,
757 : ! if surface pressures are provided (either passed as argument or in the
758 : ! HEMCO configuration file as field PSFC), pressure edges are calculated
759 : ! from PSFC and the vertical grid coordinates (Ap and Bp for a hybrid sigma
760 : ! coordinate system). The temperature field TK [K] is needed to approximate
761 : ! box heights and/or geopotential height (via the hydrostatic equation).
762 : !\\
763 : !\\
764 : ! !INTERFACE:
765 : !
766 0 : SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )
767 : !
768 : ! !USES
769 : !
770 : USE HCO_Arr_Mod, ONLY : HCO_ArrAssert, HCO_ArrCleanup
771 : USE HCO_STATE_MOD, ONLY : HCO_STATE
772 : USE HCO_CALC_MOD, ONLY : HCO_EvalFld
773 : !
774 : ! !INPUT PARAMETERS:
775 : !
776 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
777 : REAL(hp), POINTER :: PSFC(:,:) ! surface pressure (Pa)
778 : REAL(hp), POINTER :: ZSFC(:,:) ! surface geopotential height (m)
779 : REAL(hp), POINTER :: TK (:,:,:) ! air temperature (K)
780 : REAL(hp), POINTER :: BXHEIGHT(:,:,:) ! grid box height (m)
781 : REAL(hp), POINTER :: PEDGE(:,:,:) ! pressure edges (Pa)
782 : !
783 : ! !INPUT/OUTPUT PARAMETERS:
784 : !
785 : INTEGER, INTENT(INOUT) :: RC
786 : !
787 : ! !REVISION HISTORY:
788 : ! 28 Sep 2015 - C. Keller - Initial version
789 : ! See https://github.com/geoschem/hemco for complete history
790 : !EOP
791 : !------------------------------------------------------------------------------
792 : !BOC
793 : !
794 : ! !LOCAL VARIABLES:
795 : !
796 : INTEGER :: NX, NY, NZ
797 : INTEGER :: I, J, L
798 : LOGICAL :: Verb
799 : LOGICAL :: FoundPSFC
800 : LOGICAL :: FoundZSFC
801 : LOGICAL :: FoundTK
802 : LOGICAL :: FoundPEDGE
803 : LOGICAL :: FoundBXHEIGHT
804 : LOGICAL :: ERRBX, ERRZSFC
805 : REAL(hp) :: P1, P2
806 0 : REAL(hp), ALLOCATABLE, TARGET :: TmpTK(:,:,:)
807 0 : REAL(hp), POINTER :: ThisTK(:,:,:)
808 : CHARACTER(LEN=255) :: MSG
809 : CHARACTER(LEN=255) :: LOC = 'HCO_CalcVertGrid (hco_geotools_mod.F90)'
810 :
811 : LOGICAL, SAVE :: FIRST = .TRUE.
812 : LOGICAL, SAVE :: EVAL_PSFC = .TRUE.
813 : LOGICAL, SAVE :: EVAL_ZSFC = .TRUE.
814 : LOGICAL, SAVE :: EVAL_TK = .TRUE.
815 : LOGICAL, SAVE :: EVAL_PEDGE = .TRUE.
816 : LOGICAL, SAVE :: EVAL_BXHEIGHT = .TRUE.
817 :
818 : !-------------------------------
819 : ! HCO_CalcVertGrid begins here
820 : !-------------------------------
821 :
822 : ! Init
823 0 : Verb = .FALSE.
824 0 : FoundPSFC = .FALSE.
825 0 : FoundZSFC = .FALSE.
826 0 : FoundTK = .FALSE.
827 0 : FoundPEDGE = .FALSE.
828 0 : FoundBXHEIGHT = .FALSE.
829 0 : ThisTK => NULL()
830 :
831 : ! Verbose statements
832 0 : IF ( HcoState%amIRoot .AND. FIRST .AND. &
833 : HCO_IsVerb(HcoState%Config%Err) ) THEN
834 : Verb = .TRUE.
835 : ENDIF
836 : IF ( Verb ) THEN
837 0 : MSG = 'Details about vertical grid calculations (only shown on first time step):'
838 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1='-')
839 0 : MSG = '1. Input data availability: '
840 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
841 : ENDIF
842 :
843 : ! ------------------------------------------------------------------
844 : ! TK
845 : ! ------------------------------------------------------------------
846 :
847 : ! If associated, make sure that array size is correct
848 : ! and pass to HEMCO surface pressure field
849 0 : IF ( ASSOCIATED(TK) ) THEN
850 0 : NX = SIZE(TK,1)
851 0 : NY = SIZE(TK,2)
852 0 : NZ = SIZE(TK,3)
853 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
854 : NZ /= HcoState%NZ ) THEN
855 0 : WRITE(MSG,*) 'Wrong TK array size: ', NX, NY, NZ, &
856 0 : '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ
857 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
858 0 : RETURN
859 : ENDIF
860 :
861 : ! TK is not a field in grid, so don't pass
862 0 : ThisTK => TK
863 0 : FoundTK = .TRUE.
864 :
865 : ! Verbose
866 0 : IF ( Verb ) THEN
867 0 : WRITE(MSG,*) ' - Temperature field TK obtained from model interface (min,max): ', MINVAL(ThisTK), MAXVAL(ThisTK)
868 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
869 : ENDIF
870 :
871 0 : ELSEIF ( EVAL_TK ) THEN
872 0 : ALLOCATE(TmpTK(HcoState%NX,HcoState%NY,HcoState%NZ))
873 0 : CALL HCO_EvalFld( HcoState, 'TK', TmpTK, RC, FOUND=FoundTK )
874 0 : IF ( RC /= HCO_SUCCESS ) THEN
875 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
876 0 : RETURN
877 : ENDIF
878 :
879 : ! TK is sometimes listed as TMPU, so look for that too (bmy, 3/5/21)
880 0 : IF ( .not. FoundTK ) THEN
881 0 : CALL HCO_EvalFld( HcoState, 'TMPU', TmpTK, RC, FOUND=FoundTK )
882 0 : IF ( RC /= HCO_SUCCESS ) THEN
883 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
884 0 : RETURN
885 : ENDIF
886 : ENDIF
887 :
888 0 : EVAL_TK = FoundTK
889 0 : IF ( FoundTK ) ThisTK => TmpTk
890 :
891 : ! Verbose
892 0 : IF ( Verb ) THEN
893 0 : IF ( FoundTK ) THEN
894 0 : WRITE(MSG,*) ' - Temperature field TK [K] obtained from configuration file (min,max): ', MINVAL(ThisTK), MAXVAL(ThisTK)
895 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
896 : ELSE
897 0 : WRITE(MSG,*) ' - No temperature field TK found - some vertical grid calculations may not be performed...'
898 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
899 : ENDIF
900 : ENDIF
901 : ENDIF
902 :
903 : ! ------------------------------------------------------------------
904 : ! PSFC
905 : ! ------------------------------------------------------------------
906 0 : CALL HCO_ArrAssert( HcoState%Grid%PSFC, HcoState%NX, HcoState%NY, RC )
907 0 : IF ( RC /= HCO_SUCCESS ) THEN
908 0 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
909 0 : RETURN
910 : ENDIF
911 :
912 : ! If associated, make sure that array size is correct
913 : ! and pass to HEMCO surface pressure field
914 0 : IF ( ASSOCIATED(PSFC) ) THEN
915 0 : NX = SIZE(PSFC,1)
916 0 : NY = SIZE(PSFC,2)
917 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
918 0 : WRITE(MSG,*) 'Wrong PSFC array size: ', NX, NY, &
919 0 : '; should be: ', HcoState%NX, HcoState%NY
920 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
921 0 : RETURN
922 : ENDIF
923 :
924 : ! Pass to HcoState array
925 0 : HcoState%Grid%PSFC%Val = PSFC
926 0 : FoundPSFC = .TRUE.
927 :
928 : ! Verbose
929 0 : IF ( Verb ) THEN
930 0 : WRITE(MSG,*) ' - Surface pressure PSFC [Pa] obtained from model interface (min, max): ', MINVAL(HcoState%Grid%PSFC%Val), MAXVAL(HcoState%Grid%PSFC%VAL)
931 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
932 : ENDIF
933 :
934 : ! Otherwise, try to read from HEMCO configuration file
935 0 : ELSEIF ( EVAL_PSFC ) THEN
936 : CALL HCO_EvalFld( HcoState, 'PSFC', HcoState%Grid%PSFC%Val, RC, &
937 0 : FOUND=FoundPSFC )
938 0 : IF ( RC /= HCO_SUCCESS ) THEN
939 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
940 0 : RETURN
941 : ENDIF
942 :
943 : ! PSFC is sometimes listed as PS, so look for that too (bmy, 3/4/21)
944 0 : IF ( .not. FoundPSFC ) THEN
945 : CALL HCO_EvalFld( HcoState, 'PS', HcoState%Grid%PSFC%Val, RC, &
946 0 : FOUND=FoundPSFC )
947 0 : IF ( RC /= HCO_SUCCESS ) THEN
948 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
949 0 : RETURN
950 : ENDIF
951 : ENDIF
952 :
953 0 : EVAL_PSFC = FoundPSFC
954 :
955 : ! Verbose
956 0 : IF ( Verb ) THEN
957 0 : IF ( FoundPSFC ) THEN
958 0 : WRITE(MSG,*) ' - Surface pressure PSFC [Pa] obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%PSFC%Val), MAXVAL(HcoState%Grid%PSFC%VAL)
959 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
960 : ELSE
961 0 : MSG = ' - Surface pressure PSFC not found. Will attempt to calculate it.'
962 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
963 : ENDIF
964 : ENDIF
965 : ENDIF
966 :
967 : ! ------------------------------------------------------------------
968 : ! ZSFC
969 : ! ------------------------------------------------------------------
970 0 : CALL HCO_ArrAssert( HcoState%Grid%ZSFC, HcoState%NX, HcoState%NY, RC )
971 0 : IF ( RC /= HCO_SUCCESS ) THEN
972 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
973 0 : RETURN
974 : ENDIF
975 :
976 : ! If associated, make sure that array size is correct
977 : ! and pass to HEMCO surface pressure field
978 0 : IF ( ASSOCIATED(ZSFC) ) THEN
979 0 : NX = SIZE(ZSFC,1)
980 0 : NY = SIZE(ZSFC,2)
981 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
982 0 : WRITE(MSG,*) 'Wrong ZSFC array size: ', NX, NY, &
983 0 : '; should be: ', HcoState%NX, HcoState%NY
984 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
985 0 : RETURN
986 : ENDIF
987 :
988 : ! Pass to HcoState array
989 0 : HcoState%Grid%ZSFC%Val = ZSFC
990 0 : FoundZSFC = .TRUE.
991 :
992 : ! Verbose
993 0 : IF ( Verb ) THEN
994 0 : WRITE(MSG,*) ' - Surface geopotential height ZSFC [m] obtained from model interface (min, max): ', MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%VAL)
995 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
996 : ENDIF
997 :
998 : ! Otherwise, try to read from HEMCO configuration file
999 0 : ELSEIF ( EVAL_ZSFC ) THEN
1000 0 : CALL HCO_EvalFld ( HcoState, 'ZSFC', HcoState%Grid%ZSFC%Val, RC, FOUND=FoundZSFC )
1001 0 : IF ( RC /= HCO_SUCCESS ) THEN
1002 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
1003 0 : RETURN
1004 : ENDIF
1005 0 : EVAL_ZSFC = FoundZSFC
1006 :
1007 : ! Verbose
1008 0 : IF ( Verb ) THEN
1009 0 : IF ( FoundZSFC ) THEN
1010 0 : WRITE(MSG,*) ' - Surface geopotential height ZSFC [m] obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%VAL)
1011 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1012 : ELSE
1013 0 : MSG = ' - Surface geopotential height ZSFC not found. Will attempt to calculate it.'
1014 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1015 : ENDIF
1016 : ENDIF
1017 : ENDIF
1018 :
1019 : ! ------------------------------------------------------------------
1020 : ! PEDGE
1021 : ! ------------------------------------------------------------------
1022 : CALL HCO_ArrAssert( HcoState%Grid%PEDGE, HcoState%NX, &
1023 0 : HcoState%NY, HcoState%NZ+1, RC )
1024 0 : IF ( RC /= HCO_SUCCESS ) THEN
1025 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
1026 0 : RETURN
1027 : ENDIF
1028 :
1029 0 : IF ( ASSOCIATED( PEDGE ) ) THEN
1030 :
1031 0 : NX = SIZE(PEDGE,1)
1032 0 : NY = SIZE(PEDGE,2)
1033 0 : NZ = SIZE(PEDGE,3)
1034 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
1035 : NZ /= (HcoState%NZ + 1) ) THEN
1036 0 : WRITE(MSG,*) 'Wrong PEDGE array size: ', NX, NY, NZ, &
1037 0 : '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ+1
1038 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1039 0 : RETURN
1040 : ENDIF
1041 :
1042 0 : HcoState%Grid%PEDGE%Val = PEDGE
1043 0 : FoundPEDGE = .TRUE.
1044 :
1045 : ! Verbose
1046 0 : IF ( Verb ) THEN
1047 0 : WRITE(MSG,*) ' - Pressure edges PEDGE obtained from model interface (min, max): ', MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%VAL)
1048 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1049 : ENDIF
1050 :
1051 0 : ELSEIF ( EVAL_PEDGE ) THEN
1052 : CALL HCO_EvalFld ( HcoState, 'PEDGE', &
1053 0 : HcoState%Grid%PEDGE%Val, RC, FOUND=FoundPEDGE )
1054 0 : IF ( RC /= HCO_SUCCESS ) THEN
1055 0 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
1056 0 : RETURN
1057 : ENDIF
1058 0 : EVAL_PEDGE = FoundPEDGE
1059 :
1060 : ! Verbose
1061 0 : IF ( Verb ) THEN
1062 0 : IF ( FoundPEDGE ) THEN
1063 0 : WRITE(MSG,*) ' - Pressure edges PEDGE obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%VAL)
1064 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1065 : ELSE
1066 0 : MSG = ' - Pressure edges PEDGE not found. Will attempt to calculate it.'
1067 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1068 : ENDIF
1069 : ENDIF
1070 : ENDIF
1071 :
1072 : ! ------------------------------------------------------------------
1073 : ! BXHEIGHT
1074 : ! ------------------------------------------------------------------
1075 : CALL HCO_ArrAssert( HcoState%Grid%BXHEIGHT_M, HcoState%NX, &
1076 0 : HcoState%NY, HcoState%NZ, RC )
1077 0 : IF ( RC /= HCO_SUCCESS ) THEN
1078 0 : CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
1079 0 : RETURN
1080 : ENDIF
1081 :
1082 0 : IF ( ASSOCIATED( BXHEIGHT ) ) THEN
1083 :
1084 0 : NX = SIZE(BXHEIGHT,1)
1085 0 : NY = SIZE(BXHEIGHT,2)
1086 0 : NZ = SIZE(BXHEIGHT,3)
1087 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
1088 : NZ /= HcoState%NZ ) THEN
1089 0 : WRITE(MSG,*) 'Wrong BXHEIGHT array size: ', NX, NY, NZ, &
1090 0 : '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ
1091 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1092 0 : RETURN
1093 : ENDIF
1094 :
1095 0 : HcoState%Grid%BXHEIGHT_M%Val = BXHEIGHT
1096 0 : FoundBXHEIGHT = .TRUE.
1097 :
1098 : ! Verbose
1099 0 : IF ( Verb ) THEN
1100 0 : WRITE(MSG,*) ' - Boxheights BXHEIGHT_M obtained from model interface (min, max): ', MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%VAL)
1101 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1102 : ENDIF
1103 :
1104 : ! Otherwise, try to read from HEMCO configuration file
1105 0 : ELSEIF ( EVAL_BXHEIGHT ) THEN
1106 : CALL HCO_EvalFld ( HcoState, 'BXHEIGHT_M', &
1107 0 : HcoState%Grid%BXHEIGHT_M%Val, RC, FOUND=FoundBXHEIGHT )
1108 0 : IF ( RC /= HCO_SUCCESS ) THEN
1109 0 : CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
1110 0 : RETURN
1111 : ENDIF
1112 0 : EVAL_BXHEIGHT = FoundBXHEIGHT
1113 :
1114 : ! Verbose
1115 0 : IF ( Verb ) THEN
1116 0 : IF ( FoundBXHEIGHT ) THEN
1117 0 : WRITE(MSG,*) ' - Boxheights BXHEIGHT_M obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%VAL)
1118 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1119 : ELSE
1120 0 : MSG = ' - Boxheights BXHEIGHT_M not found. Will attempt to calculate it.'
1121 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1122 : ENDIF
1123 : ENDIF
1124 : ENDIF
1125 :
1126 : ! If Box height isn't available, free its memory
1127 : ! NOTE: Using NULL instead of deallocate here causes a memory leak!
1128 0 : IF ( .NOT. EVAL_BXHEIGHT .OR. .NOT. FoundBXHEIGHT ) THEN
1129 0 : CALL HCO_ArrCleanup( HcoState%Grid%BXHEIGHT_M, DeepClean=.TRUE. )
1130 : ENDIF
1131 :
1132 : ! ------------------------------------------------------------------
1133 : ! Calculate various quantities: the goal is to have the following
1134 : ! quantities defined: ZSFC, PSFC, PEDGE, BXHEIGHT_M.
1135 : ! - If PEDGE is not yet defined, it is calculated from surface
1136 : ! pressure (PSFC).
1137 : ! - If BXHEIGHT is not yet defined, it is calculated from PEDGE
1138 : ! and TK.
1139 : ! - If PSFC is not yet defined, it is set to the first level of
1140 : ! PEDGE (if defined) or uniformly set to 101325 Pa.
1141 : ! - If ZSFC is not yet defined, it is calculated from PSFC and TK.
1142 : ! - If TK is not defined, no attempt is made to initialize it to
1143 : ! a useful quantity. Calculations that require TK are omitted.
1144 : ! ------------------------------------------------------------------
1145 :
1146 : ! Verbose
1147 0 : IF ( Verb ) THEN
1148 0 : MSG = '2. Grid calculations: '
1149 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
1150 : ENDIF
1151 :
1152 : ! Set PSFC
1153 0 : IF ( .NOT. FoundPSFC ) THEN
1154 0 : IF ( FoundPEDGE ) THEN
1155 0 : HcoState%Grid%PSFC%Val(:,:) = HcoState%Grid%PEDGE%Val(:,:,1)
1156 :
1157 : ! Verbose
1158 0 : IF ( Verb ) THEN
1159 0 : MSG = ' - Surface pressure set to surface pressure edge.'
1160 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1161 : ENDIF
1162 : ELSE
1163 0 : HcoState%Grid%PSFC%Val(:,:) = 101325.0_hp
1164 0 : IF ( FIRST .AND. HcoState%amIRoot ) THEN
1165 : MSG = 'Surface pressure PSFC uniformly set to 101325 Pa! ' // &
1166 : 'This may affect the accuracy of vertical grid ' // &
1167 : 'quantities. It is recommended you provide PSFC via '// &
1168 0 : 'the model-HEMCO interface or the HEMCO configuration file!'
1169 0 : CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
1170 : ENDIF
1171 :
1172 : ! Verbose
1173 0 : IF ( Verb ) THEN
1174 0 : MSG = ' - Surface pressure uniformly set to 101325.0 Pa.'
1175 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1176 : ENDIF
1177 : ENDIF
1178 0 : FoundPSFC = .TRUE.
1179 : ENDIF
1180 :
1181 : ! Set PEDGE
1182 0 : IF ( .NOT. FoundPEDGE ) THEN
1183 : !$OMP PARALLEL DO &
1184 : !$OMP DEFAULT( SHARED ) &
1185 : !$OMP PRIVATE( I, J, L )
1186 0 : DO L = 1, HcoState%NZ+1
1187 0 : DO J = 1, HcoState%NY
1188 0 : DO I = 1, HcoState%NX
1189 0 : HcoState%Grid%PEDGE%Val(I,J,L) &
1190 0 : = HcoState%Grid%zGrid%AP(L) &
1191 0 : + ( HcoState%Grid%zGrid%BP(L) &
1192 0 : * HcoState%Grid%PSFC%Val(I,J) )
1193 : ENDDO
1194 : ENDDO
1195 : ENDDO
1196 : !$OMP END PARALLEL DO
1197 0 : FoundPEDGE = .TRUE.
1198 :
1199 : ! Verbose
1200 0 : IF ( Verb ) THEN
1201 0 : WRITE(MSG,*) ' - PEDGE calculated from PSFC, Ap, and Bp (min, max): ', &
1202 0 : MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%Val)
1203 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1204 : ENDIF
1205 : ENDIF
1206 :
1207 : ! Set surface height and/or grid box height
1208 0 : IF ( .NOT. FoundZSFC .OR. .NOT. FoundBXHEIGHT ) THEN
1209 0 : IF ( FoundTK .AND. FoundPEDGE ) THEN
1210 :
1211 : ! Initialize error flags
1212 0 : ERRZSFC = .FALSE.
1213 0 : ERRBX = .FALSE.
1214 :
1215 : ! Make sure box height is initialized if it will be calculated
1216 : CALL HCO_ArrAssert( HcoState%Grid%BXHEIGHT_M, HcoState%NX, &
1217 0 : HcoState%NY, HcoState%NZ, RC )
1218 0 : IF ( RC /= HCO_SUCCESS ) THEN
1219 0 : CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
1220 0 : RETURN
1221 : ENDIF
1222 :
1223 : !$OMP PARALLEL DO &
1224 : !$OMP DEFAULT( SHARED ) &
1225 : !$OMP PRIVATE( I, J, L, P1, P2 )
1226 0 : DO L = 1, HcoState%NZ
1227 0 : DO J = 1, HcoState%NY
1228 0 : DO I = 1, HcoState%NX
1229 :
1230 : ! BOXHEIGHT (hydrostatic equation)
1231 0 : IF ( .NOT. FoundBXHEIGHT ) THEN
1232 0 : P1 = HcoState%Grid%PEDGE%Val(I,J,L)
1233 0 : P2 = HcoState%Grid%PEDGE%Val(I,J,L+1)
1234 0 : IF ( P2 == 0.0_hp ) THEN
1235 : ERRBX = .TRUE.
1236 : ELSE
1237 0 : HcoState%Grid%BXHEIGHT_M%Val(I,J,L) = HcoState%Phys%Rdg0 &
1238 0 : * ThisTK(I,J,1) &
1239 0 : * LOG( P1 / P2 )
1240 : ENDIF
1241 : ENDIF
1242 :
1243 : ! ZSFC
1244 0 : IF ( L == 1 .AND. .NOT. FoundZSFC ) THEN
1245 0 : P1 = 101325.0_hp
1246 0 : P2 = HcoState%Grid%PEDGE%Val(I,J,1)
1247 0 : IF ( P2 == 0.0_hp ) THEN
1248 : ERRZSFC = .TRUE.
1249 : ELSE
1250 0 : HcoState%Grid%ZSFC%Val(I,J) = HcoState%Phys%Rdg0 &
1251 0 : * ThisTK(I,J,1) &
1252 0 : * LOG( P1 / P2 )
1253 : ENDIF
1254 : ENDIF
1255 : ENDDO
1256 : ENDDO
1257 : ENDDO
1258 : !$OMP END PARALLEL DO
1259 :
1260 0 : IF ( ERRZSFC ) THEN
1261 : MSG = 'Cannot calculate surface geopotential heights - at least one ' // &
1262 : 'surface pressure value is zero! You can either provide an ' // &
1263 : 'updated pressure edge field (PEDGE) or add a field with the ' // &
1264 0 : 'surface geopotential height to your configuration file (ZSFC)'
1265 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1266 0 : RETURN
1267 : ELSE
1268 0 : FoundZSFC = .TRUE.
1269 :
1270 : ! Verbose
1271 0 : IF ( Verb ) THEN
1272 0 : WRITE(MSG,*) ' - ZSFC calculated from PSFC and T (min, max): ', &
1273 0 : MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%Val)
1274 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1275 : ENDIF
1276 : ENDIF
1277 :
1278 0 : IF ( ERRBX ) THEN
1279 : MSG = 'Cannot calculate grid box heights - at least one ' // &
1280 : 'pressure value is zero! You can either provide an ' // &
1281 : 'updated pressure edge field (PEDGE) or add a field with the ' // &
1282 0 : 'grid box heights to your configuration file (BOXHEIGHT_M)'
1283 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1284 0 : RETURN
1285 : ELSE
1286 0 : FoundZSFC = .TRUE.
1287 :
1288 : ! Verbose
1289 0 : IF ( Verb ) THEN
1290 0 : WRITE(MSG,*) ' - Boxheights calculated from PEDGE and T (min, max): ', &
1291 0 : MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%Val)
1292 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1293 : ENDIF
1294 : ENDIF
1295 :
1296 : ! PEDGE and/or TK not defined
1297 : ELSE
1298 0 : IF ( .NOT. FoundZSFC .AND. FIRST .AND. HcoState%amIRoot ) THEN
1299 : MSG = 'Cannot set surface height ZSFC. This may cause ' // &
1300 : 'some extensions to fail. HEMCO tries to calculate ' // &
1301 : 'ZSFC from surface pressure and air temperature, but ' // &
1302 0 : 'at least one of these variables seem to be missing.'
1303 0 : CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
1304 : ENDIF
1305 0 : IF ( .NOT. FoundBXHEIGHT .AND. FIRST .AND. HcoState%amIRoot ) THEN
1306 : MSG = 'Cannot set boxheights BXHEIGHT_M. This may cause ' // &
1307 : 'some extensions to fail. HEMCO tries to calculate ' // &
1308 : 'BXHEIGHT from pressure edges and air temperature, but ' // &
1309 0 : 'at least one of these variables seem to be missing.'
1310 0 : CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
1311 : ENDIF
1312 : ENDIF
1313 : ENDIF
1314 :
1315 : ! ------------------------------------------------------------------
1316 : ! Wrap up and leave
1317 : ! ------------------------------------------------------------------
1318 :
1319 : ! Verbose
1320 0 : IF ( Verb ) THEN
1321 0 : WRITE(MSG,*) 'Vertical grid calculations done.'
1322 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
1323 : ENDIF
1324 :
1325 : ! Cleanup
1326 0 : ThisTK => NULL()
1327 0 : IF ( ALLOCATED(TmpTK) ) DEALLOCATE(TmpTK)
1328 :
1329 : ! Update first flag
1330 0 : FIRST = .FALSE.
1331 :
1332 : ! Return w/ success
1333 0 : RC = HCO_SUCCESS
1334 :
1335 0 : END SUBROUTINE HCO_CalcVertGrid
1336 : !EOC
1337 : !------------------------------------------------------------------------------
1338 : ! Harmonized Emissions Component (HEMCO) !
1339 : !------------------------------------------------------------------------------
1340 : !BOP
1341 : !
1342 : ! !IROUTINE: HCO_SetPBLm
1343 : !
1344 : ! !DESCRIPTION: Subroutine HCO\_SetPBLm sets the HEMCO PBL mixing height in
1345 : ! meters. It first tries to read it from field 'FldName' (from the HEMCO data
1346 : ! list), then to fill it from field 'PBLM', and then assigns the default value
1347 : ! 'DefVal' to it.
1348 : !\\
1349 : !\\
1350 : ! !INTERFACE:
1351 : !
1352 0 : SUBROUTINE HCO_SetPBLm( HcoState, FldName, PBLM, DefVal, RC )
1353 : !
1354 : ! !USES
1355 : !
1356 : USE HCO_Arr_Mod, ONLY : HCO_ArrAssert
1357 : USE HCO_STATE_MOD, ONLY : HCO_STATE
1358 : USE HCO_CALC_MOD, ONLY : HCO_EvalFld
1359 : !
1360 : ! !INPUT PARAMETERS:
1361 : !
1362 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1363 : CHARACTER(LEN=*), OPTIONAL, INTENT(IN ) :: FldName ! field name
1364 : REAL(hp), OPTIONAL, POINTER :: PBLM(:,:) ! pbl mixing height
1365 : REAL(hp), OPTIONAL, INTENT(IN ) :: DefVal ! default value
1366 : !
1367 : ! !INPUT/OUTPUT PARAMETERS:
1368 : !
1369 : INTEGER, INTENT(INOUT) :: RC
1370 : !
1371 : ! !REVISION HISTORY:
1372 : ! 28 Sep 2015 - C. Keller - Initial version
1373 : ! See https://github.com/geoschem/hemco for complete history
1374 : !EOP
1375 : !------------------------------------------------------------------------------
1376 : !BOC
1377 : !
1378 : ! !LOCAL VARIABLES:
1379 : !
1380 : INTEGER :: NX, NY
1381 : LOGICAL :: FOUND
1382 : CHARACTER(LEN=255) :: MSG
1383 : CHARACTER(LEN=255) :: LOC = 'HCO_SetPBLm (hco_geotools_mod.F90)'
1384 :
1385 : !-------------------------------
1386 : ! HCO_SetPBLm begins here
1387 : !-------------------------------
1388 :
1389 : ! Init
1390 0 : FOUND = .FALSE.
1391 :
1392 : ! Make sure size is ok. Allocate if unallocated.
1393 : CALL HCO_ArrAssert( HcoState%Grid%PBLHEIGHT, &
1394 0 : HcoState%NX, HcoState%NY, RC )
1395 0 : IF ( RC /= HCO_SUCCESS ) THEN
1396 0 : CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
1397 0 : RETURN
1398 : ENDIF
1399 :
1400 : ! Try to read from file first
1401 0 : IF ( PRESENT( FldName ) ) THEN
1402 : CALL HCO_EvalFld ( HcoState, FldName, &
1403 0 : HcoState%Grid%PBLHEIGHT%Val, RC, FOUND=FOUND )
1404 0 : IF ( RC /= HCO_SUCCESS ) THEN
1405 0 : CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
1406 0 : RETURN
1407 : ENDIF
1408 :
1409 : ! Verbose
1410 0 : IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
1411 0 : IF ( FOUND ) THEN
1412 0 : WRITE(MSG,*) 'HEMCO PBL heights obtained from field ',TRIM(FldName)
1413 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
1414 : ENDIF
1415 : ENDIF
1416 : ENDIF
1417 :
1418 : ! Pass 2D field if available
1419 0 : IF ( .not. FOUND ) THEN
1420 0 : IF ( PRESENT( PBLM ) ) THEN
1421 0 : IF ( ASSOCIATED(PBLM) ) THEN
1422 0 : NX = SIZE(PBLM,1)
1423 0 : NY = SIZE(PBLM,2)
1424 0 : IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
1425 0 : WRITE(MSG,*) 'Wrong PBLM array size: ', NX, NY, &
1426 0 : '; should be: ', HcoState%NX, HcoState%NY
1427 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1428 0 : RETURN
1429 : ENDIF
1430 :
1431 : ! Pass data
1432 0 : HcoState%Grid%PBLHEIGHT%Val = PBLM
1433 0 : FOUND = .TRUE.
1434 :
1435 : ! Verbose
1436 0 : IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
1437 0 : WRITE(MSG,*) 'HEMCO PBL heights obtained from provided 2D field.'
1438 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
1439 : ENDIF
1440 : ENDIF
1441 : ENDIF
1442 : ENDIF
1443 :
1444 : ! Finally, assign default value if field not yet set
1445 : ! IF ( .NOT. FOUND .AND. PRESENT(DefVal) ) THEN
1446 0 : IF ( .NOT. FOUND ) THEN
1447 0 : IF ( PRESENT(DefVal) ) THEN
1448 : ! Make sure size is ok
1449 : CALL HCO_ArrAssert( HcoState%Grid%PBLHEIGHT, &
1450 0 : HcoState%NX, HcoState%NY, RC )
1451 0 : IF ( RC /= HCO_SUCCESS ) THEN
1452 0 : CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
1453 0 : RETURN
1454 : ENDIF
1455 :
1456 : ! Pass data
1457 0 : HcoState%Grid%PBLHEIGHT%Val = DefVal
1458 0 : FOUND = .TRUE.
1459 :
1460 : ! Verbose
1461 0 : IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
1462 0 : WRITE(MSG,*) 'HEMCO PBL heights uniformly set to ', DefVal
1463 0 : CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
1464 : ENDIF
1465 : ENDIF
1466 : ENDIF
1467 :
1468 : ! Error check
1469 0 : IF ( .NOT. FOUND ) THEN
1470 0 : WRITE(MSG,*) 'Cannot set PBL height: a valid HEMCO data field, ', &
1471 0 : 'an explicit 2D field or a default value must be provided!'
1472 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1473 0 : RETURN
1474 : ENDIF
1475 :
1476 : ! Return w/ success
1477 0 : RC = HCO_SUCCESS
1478 :
1479 : END SUBROUTINE HCO_SetPBLm
1480 : !EOC
1481 : !!------------------------------------------------------------------------------
1482 : !! Harmonized Emissions Component (HEMCO) !
1483 : !!------------------------------------------------------------------------------
1484 : !!BOP
1485 : !!
1486 : !! !SUBROUTINE: HCO_CalcPBLlev3D
1487 : !!
1488 : !! !DESCRIPTION:
1489 : !!\\
1490 : !!\\
1491 : !! !INTERFACE:
1492 : !!
1493 : ! SUBROUTINE HCO_CalcPBLlev3D ( HcoState, PBLFRAC, RC )
1494 : !!
1495 : !! !USES
1496 : !!
1497 : ! USE HCO_Arr_Mod, ONLY : HCO_ArrAssert
1498 : ! USE HCO_STATE_MOD, ONLY : HCO_STATE
1499 : !!
1500 : !! !INPUT PARAMETERS:
1501 : !!
1502 : ! TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1503 : ! REAL(hp), POINTER :: PBLFRAC(:,:,:) ! planetary PBL fraction
1504 : !!
1505 : !! !INPUT/OUTPUT PARAMETERS:
1506 : !!
1507 : ! INTEGER, INTENT(INOUT) :: RC
1508 : !!
1509 : !! !REVISION HISTORY:
1510 : !! 05 May 2016 - C. Keller - Initial version
1511 : !! See https://github.com/geoschem/hemco for complete history
1512 : !!EOP
1513 : !!------------------------------------------------------------------------------
1514 : !!BOC
1515 : !!
1516 : !! !LOCAL VARIABLES:
1517 : !!
1518 : ! INTEGER :: I, J, L
1519 : ! CHARACTER(LEN=255) :: MSG
1520 : ! CHARACTER(LEN=255) :: LOC = 'HCO_CalcPBLlev3D (hco_geotools_mod.F90)'
1521 : !
1522 : ! !-------------------------------
1523 : ! ! HCO_CalcPBLlev3D begins here
1524 : ! !-------------------------------
1525 : !
1526 : ! ! Check input array size
1527 : ! IF ( SIZE(PBLFRAC,1) /= HcoState%NX .OR. &
1528 : ! SIZE(PBLFRAC,2) /= HcoState%NY .OR. &
1529 : ! SIZE(PBLFRAC,3) /= HcoState%NZ ) THEN
1530 : ! WRITE(MSG,*) 'Input array PBLFRAC has wrong horiz. dimensions: ', &
1531 : ! SIZE(PBLFRAC,1),SIZE(PBLFRAC,2),SIZE(PBLFRAC,3)
1532 : ! CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1533 : ! RETURN
1534 : ! ENDIF
1535 : !
1536 : ! ! Make sure array is associated
1537 : ! CALL HCO_ArrAssert( HcoState%Grid%PBL, HcoState%NX, HcoState%NY, RC )
1538 : ! IF ( RC /= HCO_SUCCESS ) RETURN
1539 : !
1540 : ! ! Initialize values
1541 : ! HcoState%Grid%PBL%Val = 1.0
1542 : !
1543 : !!$OMP PARALLEL DO &
1544 : !!$OMP DEFAULT( SHARED ) &
1545 : !!$OMP PRIVATE( I, J, L ) &
1546 : !!$OMP SCHEDULE( DYNAMIC )
1547 : ! DO J = 1, HcoState%NY
1548 : ! DO I = 1, HcoState%NX
1549 : ! ! Search for first level where PBL fraction is zero
1550 : ! DO L = 1, HcoState%NZ
1551 : ! IF ( PBLFRAC(I,J,L) > 0.0_hp ) CYCLE
1552 : ! HcoState%Grid%PBL%Val(I,J) = MAX(L-1,1)
1553 : ! EXIT
1554 : ! ENDDO
1555 : ! ENDDO
1556 : ! ENDDO
1557 : !!$OMP END PARALLEL DO
1558 : !
1559 : ! ! Return w/ success
1560 : ! RC = HCO_SUCCESS
1561 : !
1562 : ! END SUBROUTINE HCO_CalcPBLlev3D
1563 : !!EOC
1564 : !!------------------------------------------------------------------------------
1565 : !! Harmonized Emissions Component (HEMCO) !
1566 : !!------------------------------------------------------------------------------
1567 : !!BOP
1568 : !!
1569 : !! !SUBROUTINE: HCO_CalcPBLlev2D
1570 : !!
1571 : !! !DESCRIPTION:
1572 : !!\\
1573 : !!\\
1574 : !! !INTERFACE:
1575 : !!
1576 : ! SUBROUTINE HCO_CalcPBLlev2D ( HcoState, PBLlev, RC )
1577 : !!
1578 : !! !USES
1579 : !!
1580 : ! USE HCO_Arr_Mod, ONLY : HCO_ArrAssert
1581 : ! USE HCO_STATE_MOD, ONLY : HCO_STATE
1582 : !!
1583 : !! !INPUT PARAMETERS:
1584 : !!
1585 : ! TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
1586 : ! INTEGER, POINTER :: PBLlev(:,:) ! planetary PBL lev
1587 : !!
1588 : !! !INPUT/OUTPUT PARAMETERS:
1589 : !!
1590 : ! INTEGER, INTENT(INOUT) :: RC
1591 : !!
1592 : !! !REVISION HISTORY:
1593 : !! 05 May 2016 - C. Keller - Initial version
1594 : !! See https://github.com/geoschem/hemco for complete history
1595 : !!EOP
1596 : !!------------------------------------------------------------------------------
1597 : !!BOC
1598 : !!
1599 : !! !LOCAL VARIABLES:
1600 : !!
1601 : ! INTEGER :: I, J, L
1602 : ! CHARACTER(LEN=255) :: MSG
1603 : ! CHARACTER(LEN=255) :: LOC = 'HCO_CalcPBLlev2D (hco_geotools_mod.F90)'
1604 : !
1605 : ! !-------------------------------
1606 : ! ! HCO_CalcPBLlev2D begins here
1607 : ! !-------------------------------
1608 : !
1609 : ! ! Make sure array is associated
1610 : ! CALL HCO_ArrAssert( HcoState%Grid%PBL, HcoState%NX, HcoState%NY, RC )
1611 : ! IF ( RC /= HCO_SUCCESS ) RETURN
1612 : !
1613 : ! ! Check input array size
1614 : ! IF ( SIZE(PBLlev,1) /= HcoState%NX .OR. SIZE(PBLlev,2) /= HcoState%NY ) THEN
1615 : ! WRITE(MSG,*) 'Input array PBLlev has wrong horiz. dimensions: ', &
1616 : ! SIZE(PBLlev,1),SIZE(PBLlev,2)
1617 : ! CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1618 : ! RETURN
1619 : ! ENDIF
1620 : !
1621 : ! ! Set values from PBLlev
1622 : !!$OMP PARALLEL DO &
1623 : !!$OMP DEFAULT( SHARED ) &
1624 : !!$OMP PRIVATE( I, J ) &
1625 : !!$OMP SCHEDULE( DYNAMIC )
1626 : ! DO J = 1, HcoState%NY
1627 : ! DO I = 1, HcoState%NX
1628 : ! HcoState%Grid%PBL%Val(I,J) = MAX(PBLlev(I,J),1)
1629 : ! ENDDO
1630 : ! ENDDO
1631 : !!$OMP END PARALLEL DO
1632 : !
1633 : ! ! Return w/ success
1634 : ! RC = HCO_SUCCESS
1635 : !
1636 : ! END SUBROUTINE HCO_CalcPBLlev2D
1637 : !!EOC
1638 : END MODULE HCO_GeoTools_Mod
1639 : !EOM
|
