Calculate huss and other variables on the fly. COSIMA/access-om2#242

COSIMA · Sep 14, 2021 · 7f2a4a0 · aekiss · Dec 2, 2021 · 7f2a4a0
1 parent f309fc9
commit 7f2a4a0
Show file tree

Hide file tree

Showing 3 changed files with 69 additions and 11 deletions.
diff --git a/libforcing/src/forcing_field.F90 b/libforcing/src/forcing_field.F90
@@ -172,18 +172,66 @@ subroutine forcing_field_calculate(self, file_index, result_array)
     integer, intent(in) :: file_index
     real, dimension(:, :), intent(inout) :: result_array
 
+    real, dimension(:, :), allocatable, target :: tmp1, tmp2
+    real, dimension(:, :), allocatable :: E
+    real, dimension(:, :), pointer :: Td, sp
+
+    integer :: nx, ny
+    real, parameter :: RDRY = 287.0597
+    real, parameter :: RVAP = 461.5250
+    real, parameter :: A1 = 611.21
+    real, parameter :: A3 = 17.502
+    real, parameter :: A4 = 32.19
+    real, parameter :: T0 = 273.16
+
     if (trim(self%product_name) == 'ERA5') then
+
+        nx = self%ncvars(1)%nx
+        ny = self%ncvars(1)%ny
+        allocate(tmp1(nx, ny), tmp2(nx, ny))
+
         if (trim(self%coupling_name) == 'rain_ai') then
-            ! Rain is calculated as total precipitation - snow
-            ! FIXME: do calculation with field name checks
-            call self%ncvars(1)%read_data(file_index, result_array)
+            ! Rain is calculated as mcpr
+            ! (mean convective precipitation rate [kg m**-2 s**-1]) plus
+            ! mlspr (mean large-scale precipitation rate [kg m**-2 s**-1])
+            call self%ncvars(1)%read_data(file_index, tmp1)
+            call self%ncvars(2)%read_data(file_index, tmp2)
+            result_array(:, :) = tmp1(:, :) + tmp2(:, :)
+
+        elseif (trim(self%coupling_name) == 'runof_ai') then
+            ! Runoff is calculated as msror
+            ! (mean surface runoff rate [kg m**-2 s**-1]) plus
+            ! mssror (mean sub-surface runoff rate [kg m**-2 s**-1])
+
+            call self%ncvars(1)%read_data(file_index, tmp1)
+            call self%ncvars(2)%read_data(file_index, tmp2)
+            result_array(:, :) = tmp1(:, :) + tmp2(:, :)
+
         elseif (trim(self%coupling_name) == 'qair_ai') then
-            ! Humidity is calculated using surface temperature and dew point
-            ! FIXME: do calculation with field name checks
-            call self%ncvars(1)%read_data(file_index, result_array)
+            ! Specific humidity at 2m
+
+            call assert(self%ncvars(1)%name == 'd2m', &
+                        'Unexpected name for in huss calculation')
+            call assert(self%ncvars(2)%name == 'sp', &
+                        'Unexpected name for surface pressure')
+
+            call self%ncvars(1)%read_data(file_index, tmp1)
+            call self%ncvars(2)%read_data(file_index, tmp2)
+            Td => tmp1
+            sp => tmp2
+
+            ! Saturation water vapour from Teten's formula
+            allocate(E(nx, ny))
+            E = A1*exp(A3*(Td-T0)/(Td-A4))
+
+            ! Saturation specific humidity at 2m
+            result_array(:, :) = (RDRY/RVAP)*E / (sp-((1-RDRY/RVAP)*E))
+            deallocate(E)
         else
             call self%ncvars(1)%read_data(file_index, result_array)
         endif
+
+        deallocate(tmp1, tmp2)
     elseif (trim(self%product_name) == 'JRA55-do') then
         call self%ncvars(1)%read_data(file_index, result_array)
     else

diff --git a/libforcing/src/util.F90 b/libforcing/src/util.F90
@@ -51,7 +51,8 @@ subroutine read_data(ncid, varid, varname, indx, dataout)
 
     integer, dimension(:), allocatable :: count, start
     real, dimension(1) :: scalar_dataout
-    integer :: ndims, nx, ny, time
+    integer :: ndims, nx, ny, time, status
+    real :: scale_factor, offset
 
     call get_var_dims(ncid, varid, ndims, nx, ny, time)
     call assert(ndims == 1 .or. ndims == 2 .or. ndims == 3 .or. ndims == 4, &
@@ -85,6 +86,15 @@ subroutine read_data(ncid, varid, varname, indx, dataout)
                     'Get var '//trim(varname))
     endif
 
+    status = nf90_get_att(ncid, varid, "scale_factor", scale_factor)
+    if (status == nf90_noerr) then
+        dataout(:, :) = dataout(:, :) * scale_factor
+    endif
+    status = nf90_get_att(ncid, varid, "add_offset", offset)
+    if (status == nf90_noerr) then
+        dataout(:, :) = dataout(:, :) + offset
+    endif
+
     deallocate(count, start)
 
 end subroutine read_data

diff --git a/tests/ERA5/forcing.json b/tests/ERA5/forcing.json
@@ -82,13 +82,13 @@
     {
         "coupling_field_name": "qair_ai",
         "input_fields": [
-            {
-                "filename": "/g/data/rt52/era5/single-levels/reanalysis/2t/{{year}}/2t_era5_oper_sfc_{{year}}{{month}}{{start_day}}-{{year}}{{month}}{{end_day}}.nc",
-                "fieldname": "t2m"
-            },
             {
                 "filename": "/g/data/rt52/era5/single-levels/reanalysis/2d/{{year}}/2d_era5_oper_sfc_{{year}}{{month}}{{start_day}}-{{year}}{{month}}{{end_day}}.nc",
                 "fieldname": "d2m"
+            },
+            {
+                "filename": "/g/data/rt52/era5/single-levels/reanalysis/sp/{{year}}/sp_era5_oper_sfc_{{year}}{{month}}{{start_day}}-{{year}}{{month}}{{end_day}}.nc",
+                "fieldname": "sp"
             }
         ]
     },