| | 1 | == Deployment of Vtables for NCEP Inputs == |
| | 2 | |
| | 3 | As with the ECMWF Vtable work of Work Order 06, the motivation for using Vtables was to facilitate a structured approach to make modifications to the GRIB messages included or excluded in a GRIB file used as input to FLEXPART. We needed a flexible approach to quickly change GRIB codes used within the source code. Prior to our work, ''readwind_gfs()'' had a general structure in which all of the GRIB codes were explicitly written in the routine and, further, there were confusing transformations done from any GRIB2 codes to GRIB1 codes, which the rest of the routine was coded for. The original FPv9.2 GRIB2 code looked as follows: |
| | 4 | |
| | 5 | |
| | 6 | {{{ |
| | 7 | |
| | 8 | !read the grib2 identifiers |
| | 9 | call grib_get_int(igrib,'discipline',discipl,iret) |
| | 10 | call grib_check(iret,gribFunction,gribErrorMsg) |
| | 11 | call grib_get_int(igrib,'parameterCategory',parCat,iret) |
| | 12 | call grib_check(iret,gribFunction,gribErrorMsg) |
| | 13 | call grib_get_int(igrib,'parameterNumber',parNum,iret) |
| | 14 | call grib_check(iret,gribFunction,gribErrorMsg) |
| | 15 | call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret) |
| | 16 | call grib_check(iret,gribFunction,gribErrorMsg) |
| | 17 | call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', & |
| | 18 | valSurf,iret) |
| | 19 | call grib_check(iret,gribFunction,gribErrorMsg) |
| | 20 | |
| | 21 | |
| | 22 | !convert to grib1 identifiers |
| | 23 | isec1(6)=-1 |
| | 24 | isec1(7)=-1 |
| | 25 | isec1(8)=-1 |
| | 26 | if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.100)) then ! T |
| | 27 | isec1(6)=11 ! indicatorOfParameter |
| | 28 | isec1(7)=100 ! indicatorOfTypeOfLevel |
| | 29 | isec1(8)=valSurf/100 ! level, convert to hPa |
| | 30 | elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U |
| | 31 | isec1(6)=33 ! indicatorOfParameter |
| | 32 | isec1(7)=100 ! indicatorOfTypeOfLevel |
| | 33 | isec1(8)=valSurf/100 ! level, convert to hPa |
| | 34 | . |
| | 35 | . |
| | 36 | . |
| | 37 | }}} |
| | 38 | |
| | 39 | |
| | 40 | and then, further down, the GRIB1 codes were checked for storing the message data: |
| | 41 | |
| | 42 | |
| | 43 | {{{ |
| | 44 | if (isec1(6).ne.-1) then |
| | 45 | |
| | 46 | do j=0,nymin1 |
| | 47 | do i=0,nxfield-1 |
| | 48 | if((isec1(6).eq.011).and.(isec1(7).eq.100)) then |
| | 49 | ! TEMPERATURE |
| | 50 | if((i.eq.0).and.(j.eq.0)) then |
| | 51 | do ii=1,nuvz |
| | 52 | if ((isec1(8)*100.0).eq.akz(ii)) numpt=ii |
| | 53 | !print *, 'isec1(8), akz(ii), numpt: ', isec1(8), akz(ii), numpt |
| | 54 | end do |
| | 55 | !stop |
| | 56 | endif |
| | 57 | help=zsec4(nxfield*(ny-j-1)+i+1) |
| | 58 | if(i.le.i180) then |
| | 59 | tth(i179+i,j,numpt,n)=help |
| | 60 | else |
| | 61 | tth(i-i181,j,numpt,n)=help |
| | 62 | endif |
| | 63 | endif |
| | 64 | if((isec1(6).eq.033).and.(isec1(7).eq.100)) then |
| | 65 | ! U VELOCITY |
| | 66 | if((i.eq.0).and.(j.eq.0)) then |
| | 67 | do ii=1,nuvz |
| | 68 | if ((isec1(8)*100.0).eq.akz(ii)) numpu=ii |
| | 69 | end do |
| | 70 | endif |
| | 71 | help=zsec4(nxfield*(ny-j-1)+i+1) |
| | 72 | if(i.le.i180) then |
| | 73 | uuh(i179+i,j,numpu)=help |
| | 74 | else |
| | 75 | uuh(i-i181,j,numpu)=help |
| | 76 | endif |
| | 77 | endif |
| | 78 | . |
| | 79 | . |
| | 80 | . |
| | 81 | }}} |
| | 82 | |
| | 83 | |
| | 84 | Clearly, a lot of numbers are hard-coded in the above, making it very precarious to consider current and future modification needs. By using the concept of a Variable Table (or Vtable) from WRF, we were able to abstract away all of the hard-coded numbers, clean up the code substantially, and make it more structured for future modifications. |
| | 85 | |
| | 86 | Before moving ahead, please note the inefficient structure of the ''i,j'' loops in the above code. At every iteration of this doubly-nested loop, a series of if-statements is accounted before actually assigning a single number to a single array element. A much more efficient approach would be to have the doubly-nested loops coded for each variable, separately, so that they are taking better advantage of cache and working more on the simple reading of array data. This is described a bit more below in our new implementation. |
| | 87 | |
| | 88 | The Vtable for NCEP GRIB2 files looks as follows |
| | 89 | |
| | 90 | |
| | 91 | {{{ |
| | 92 | fpname | paramId | indofParam | discipline | category | number | typeSurf | typeofLevel | units | shortName | description | gribVersion | center | |
| | 93 | ---------+---------+------------+------------+----------+--------+-----------+--------------------+-----------+-----------+--------------------+------------+---------+ |
| | 94 | TT | | | 0 | 0 | 0 | 100 | isobaricInhPa | k | t | temperature | 2 | kwbc | |
| | 95 | UU | | | 0 | 2 | 2 | 100 | isobaricInhPa | ms**-1 | u | u wind | 2 | kwbc | |
| | 96 | VV | | | 0 | 2 | 3 | 100 | isobaricInhPa | ms**-1 | v | v wind | 2 | kwbc | |
| | 97 | RH | | | 0 | 1 | 1 | 100 | isobaricInhPa | % | r | relative humidity | 2 | kwbc | |
| | 98 | PS | | | 0 | 3 | 0 | 1 | surface | Pa | sp | surface pressure | 2 | kwbc | |
| | 99 | WW | | | 0 | 2 | 8 | 100 | isobaricInhPa | Pa s**-1 | w | vertical velocity | 2 | kwbc | |
| | 100 | SD | | | 0 | 1 | 13 | 1 | surface | kg m-2 | sdwe | snow depth water eq| 2 | kwbc | |
| | 101 | SLP | | | 0 | 3 | 1 | 101 | meanSea | Pa | prmsl | mean sea lev. pres.| 2 | kwbc | |
| | 102 | U10 | | | 0 | 2 | 2 | 103 | heightAboveGround | ms**-1 | 10u | 10 m u wind | 2 | kwbc | |
| | 103 | V10 | | | 0 | 2 | 3 | 103 | heightAboveGround | ms**-1 | 10u | 10 m u wind | 2 | kwbc | |
| | 104 | T2 | | | 0 | 0 | 0 | 103 | heightAboveGround | k | 2t | 2m temperature | 2 | kwbc | |
| | 105 | TD2 | | | 0 | 0 | 6 | 103 | heightAboveGround | k | 2d | 2m dew point | 2 | kwbc | |
| | 106 | TCC | | | 0 | 6 | 1 | 244 | Low cloud layer | % | tcc | total cloud cover | 2 | kwbc | |
| | 107 | LSPREC | | | 0 | 1 | 7 | 1 | surface | kg m**-2 | tp | total precip. | 2 | kwbc | |
| | 108 | CONVPREC | | | 0 | 1 | 196 | 1 | surface | kg m**-2 | acpcp | convective precip. | 2 | kwbc | |
| | 109 | ORO | | | 0 | 3 | 5 | 1 | surface | gpm | orog | orography/geopot | 2 | kwbc | |
| | 110 | LSM | | | 2 | 0 | 0 | 1 | surface | none | lsm | land-sea mask | 2 | kwbc | |
| | 111 | HMIX | | | 0 | 3 | 196 | 1 | surface | m | hpbl | mixing height | 2 | kwbc | |
| | 112 | RH2 | | | 0 | 1 | 1 | 103 | heightAboveGround | % | r | relative humidity | 2 | kwbc | |
| | 113 | TSIG1 | | | 0 | 0 | 0 | 104 | sigma | k | t | temperature | 2 | kwbc | |
| | 114 | USIG1 | | | 0 | 2 | 2 | 104 | sigma | ms**-1 | u | u wind | 2 | kwbc | |
| | 115 | VSIG1 | | | 0 | 2 | 3 | 104 | sigma | ms**-1 | v | v wind | 2 | kwbc | |
| | 116 | ---------+---------+------------+------------+----------+--------+-----------+--------------------+-----------+-----------+--------------------+------------+--------+ |
| | 117 | # |
| | 118 | # This table constitutes the parameters that work for NCEP GRIB2 files available in the CTBTO operational |
| | 119 | # data, as well as the NCEP FNL Operational Model Global Tropospheric Analyses, continuing from |
| | 120 | # July 1999, available from http://rda.ucar.edu/datasets/ds083.2/ |
| | 121 | # |
| | 122 | # For the most part, parameters were chosen by looking through the codes used in readwind_gfs.f90, |
| | 123 | # but cross-checking was done using some of the information available at the following sites: |
| | 124 | # |
| | 125 | # * Fixed surface types and units -- http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-5.shtml |
| | 126 | # * GRIB2 fields -- http://rda.ucar.edu/datasets/ds330.0/metadata/grib2.html |
| | 127 | # * GRIB2, parameters for Discipline 2, Category 3 -- http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-2-0-3.shtml |
| | 128 | # |
| | 129 | # Note that there is a bug in the FLEPXART v9.02 readwind_gfs.f90 which results in TCC not being |
| | 130 | # read in correctly and, therefore, not affecting wet deposition. This was reported (with a |
| | 131 | # suggested fix) in Ticket #143 at flexpart.eu Trac ticketing system. |
| | 132 | # |
| | 133 | |
| | 134 | }}} |
| | 135 | |
| | 136 | |
| | 137 | The table maps FLEXPART names of variables (''fpname'') to all of the relevant information in the GRIB files so that FLEXPART source code only needs to refer to the ''fpname'' at all times. The revised ''readwind_gfs()'' now looks as follows |
| | 138 | |
| | 139 | |
| | 140 | {{{ |
| | 141 | ! Get the fpname |
| | 142 | fpname = vtable_get_fpname(igrib, my_vtable) |
| | 143 | . |
| | 144 | . |
| | 145 | . |
| | 146 | ! TEMPERATURE |
| | 147 | if( trim(fpname) .eq. 'TT' ) then |
| | 148 | do j=0,nymin1 |
| | 149 | do i=0,nxfield-1 |
| | 150 | if((i.eq.0).and.(j.eq.0)) then |
| | 151 | do ii=1,nuvz |
| | 152 | if (current_grib_level .eq. akz(ii)) numpt=ii |
| | 153 | end do |
| | 154 | endif |
| | 155 | help=zsec4(nxfield*(ny-j-1)+i+1) |
| | 156 | if(i.le.i180) then |
| | 157 | tth(i179+i,j,numpt,n)=help |
| | 158 | else |
| | 159 | tth(i-i181,j,numpt,n)=help |
| | 160 | endif |
| | 161 | end do |
| | 162 | end do |
| | 163 | endif |
| | 164 | |
| | 165 | |
| | 166 | ! U VELOCITY |
| | 167 | if( trim(fpname) .eq. 'UU' ) then |
| | 168 | do j=0,nymin1 |
| | 169 | do i=0,nxfield-1 |
| | 170 | if((i.eq.0).and.(j.eq.0)) then |
| | 171 | do ii=1,nuvz |
| | 172 | if (current_grib_level .eq. akz(ii)) numpu=ii |
| | 173 | end do |
| | 174 | endif |
| | 175 | help=zsec4(nxfield*(ny-j-1)+i+1) |
| | 176 | if(i.le.i180) then |
| | 177 | uuh(i179+i,j,numpu)=help |
| | 178 | else |
| | 179 | uuh(i-i181,j,numpu)=help |
| | 180 | endif |
| | 181 | end do |
| | 182 | end do |
| | 183 | endif |
| | 184 | . |
| | 185 | . |
| | 186 | . |
| | 187 | }}} |
| | 188 | |
| | 189 | Note that the only reference to the GRIB variables is now by ''fpname''. If, at a later date we want to use different GRIB codes, we only have to modify the Vtable, and there is no reason to modify source code. Also note that each variable now has its own i,j loop. This results in performance improvements of approximately 25%, relative to the original Vtable code which didn't have these optimizations (discussed a little more, below) |
| | 190 | |
| | 191 | Note that similar changes were also made in ''gridcheck_gfs()''. |
| | 192 | |
| | 193 | Throughout this development, model outputs were frequently checked in the testing environment to quickly catch a number of minor problems so that we could rapidly fix them before moving on. |
| | 194 | |
| | 195 | The GRIB1 Vtables were then implemented and the above code, after, again, encountering a number of minor problems, was able to work correctly in the test environment. We found, however, that there were no TCC, LSPREC, or CONVPREC fields available in any of the NCEP GRIB1 files. The files we obtained from CTBTO were bitwise identical with those obtained from |
| | 196 | |
| | 197 | == Optimization of the code == |
| | 198 | |
| | 199 | Since the Work Order 6 deployment of Vtables for ECMWF inputs, the performance had declined by approximately 25%, and this was worrisome. After initial deployment of the NCEP Vtables (without yet restructuring the loops as described above), we set up two side-by-side runs - FPv9.02 and our new code, FPvWO8 - and profiled them, and found that readwind_gfs() was taking about twice as long in FPvWO8. This makes sense, because the comparisons performed at EVERY i,j iteration were now being done on strings rather than simple numbers, adding significant overhead. By restructuring the code as described above, these string comparisons were only performed each time we read in a new GRIB message from the met input file. After performing this restructuring, we found execution time between the two versions was similar. |
| | 200 | |
| | 201 | Further, after optimizing the NCEP code we went back and modified both ''readwind()'' and ''readwind_nests()'' for ECMWF with similar improvements. |
| | 202 | |
| | 203 | There are other modifications we could make that would offer even greater efficiency, but the improved performance would likely be somewhat negligible. |
| | 204 | |
| | 205 | == Enhanced testing of the code == |
| | 206 | |
| | 207 | So far, we had been testing all of our code against the test environment, but the test environment is still immature in its ability to test wet deposition. We felt like it was prudent to run several FPv9.02 versus FPvWO8 tests against 48-hour simulations driven by both GRIB1 and GRIB2. We used tools in the testing environment to detect that there were no differences between the outputs |
| | 208 | |
| | 209 | However, we then decided we should try a 48-hour wet deposition case with NCEP GRIB2 inputs, and this led to several problems. |
| | 210 | * First, our testing environment had been coded to handle deposition fields that might contain multiple species, so it assumed 3D arrays (X, Y, Species). Unfortunately, unknown to us at the time, the pflexible routines for extracting FP output chose to give us only a 2D deposition field if there was only one species, but a 3D deposition field if there was more than one species. This took a while to figure out, but we were able to adjust our testing environment for that |
| | 211 | * A more serious problem emerged when we finally compared the outputs of FPv9.02 and FPvWO8 - they were different! After resolving a couple of errors in the Vtables, they continued to be different, and after a lot of investigation we found that the original FPv9.02 did not correctly handle the total cloud cover TCC field. In fact, it never read it, so it had zero effect on wet deposition. Meanwhile, our FPvWO8 had been coded to handle it properly, so results were different. Ultimately, we were able to isolate the problem and suggest a simple solution in FPv9.02 and open Ticket 143 in the main flexpart.eu ticketing system |
| | 212 | * The entire set of FPv9.02 and FPvWO8 test cases has been tarred, zipped (2.2 GBytes) |
| | 213 | |
| | 214 | == How to use Vtables with FLEXPART == |
| | 215 | |
| | 216 | The current version of our code, FPWO8, now uses Vtables. Although we could structure the code so that it used the old method if desired, or the new, Vtable method if desired, we currently can see no reason for doing this - the Vtable code has been tested as extensively as possible, and if future problems arise, we think it would be just as easy to fix in the Vtable code as in the older code. Coding both methods would make it somewhat uglier, and would simply result in more code to maintain. |
| | 217 | |
| | 218 | The use of Vtables is very simple from a user perspective - a collection of the following Vtables is available in the Vtables subdirectory of the source code distribution: |
| | 219 | |
| | 220 | |
| | 221 | {{{ |
| | 222 | Vtable_ecmwf_grib1 |
| | 223 | Vtable_ecmwf_grib1_2 |
| | 224 | Vtable_ecmwf_grib2 |
| | 225 | Vtable_ncep_grib1 |
| | 226 | Vtable_ncep_grib2 |
| | 227 | }}} |
| | 228 | |
| | 229 | FLEXPART will automatically look at the GRIB file and determine which of the above categories it falls into, and then it will expect to find the appropriate Vtable in the run directory (the same directory that contains ''pathnames''). It's that simple. |
| | 230 | |
| | 231 | As we move to FP10, we will need to explore whether we want to change how Vtables are used, particularly when we consider situations where we might use different Vtables for nests of different GRIB types (as far out as that seems). |
