iow_esm
/
components.flux_calculator
4
0
Fork 0
Code Issues 1 Pull Requests Projects Releases 1 Wiki Activity
Browse Source

Initial commit. 

Obtained as zip archive from Hagen Radtke on May 4th 2021
distribute-radiation-on-surface-types
Sven Karsten 2 years ago
commit
d6912bb54e
24 changed files with 3525 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 16
      build.sh
  2. 16
      build_library.sh
  3. 15
      build_local.sh
  4. 87
      src/.vscode/tasks.json
  5. 52
      src/decomp_def.F90
  6. 1241
      src/flux_calculator.F90
  7. 560
      src/flux_calculator_basic.F90
  8. 272
      src/flux_calculator_calculate.F90
  9. 169
      src/flux_calculator_io.F90
  10. 255
      src/flux_calculator_prepare.F90
  11. 72
      src/flux_lib/auxiliaries/flux_aux_vapor.F90
  12. 36
      src/flux_lib/constants/flux_constants.F90
  13. 20
      src/flux_lib/ferret.jnl
  14. 43
      src/flux_lib/flux_library.F90
  15. 71
      src/flux_lib/heat/flux_heat_latent.F90
  16. 101
      src/flux_lib/heat/flux_heat_sensible.F90
  17. BIN
      src/flux_lib/mass/.flux_mass_evap.F90.swp
  18. 85
      src/flux_lib/mass/flux_mass_evap.F90
  19. 76
      src/flux_lib/momentum/flux_momentum.F90
  20. 46
      src/flux_lib/radiation/flux_radiation_blackbody.F90
  21. 61
      src/function_sent.F90
  22. 112
      src/read_dimgrid.F90
  23. 98
      src/read_grid.F90
  24. 21
      src/routine_hdlerr.F90

16
build.sh
Unescape View File

@ -0,0 +1,16 @@
FC=mpiifort
mkdir build
cd build
$FC -c ../src/flux_lib/constants/*.F90
$FC -c ../src/flux_lib/auxiliaries/*.F90
$FC -c ../src/flux_lib/mass/*.F90
$FC -c ../src/flux_lib/heat/*.F90
$FC -c ../src/flux_lib/momentum/*.F90
$FC -c ../src/flux_lib/*.F90
#$FC -shared -o flux_library.so ../src/flux_lib/*.F90
cd ..

16
build_library.sh
Unescape View File

@ -0,0 +1,16 @@
FC=mpiifort
mkdir build
cd build
$FC -c ../src/flux_lib/constants/*.F90
$FC -c ../src/flux_lib/auxiliaries/*.F90
$FC -c ../src/flux_lib/mass/*.F90
$FC -c ../src/flux_lib/heat/*.F90
$FC -c ../src/flux_lib/momentum/*.F90
$FC -c ../src/flux_lib/*.F90
#$FC -shared -o flux_library.so ../src/flux_lib/*.F90
cd ..

15
build_local.sh
Unescape View File

@ -0,0 +1,15 @@
FC=gfortran
rm -rf build_local
mkdir build_local
cd build_local
$FC -c ../src/flux_lib/constants/*.F90
$FC -c ../src/flux_lib/auxiliaries/*.F90
$FC -c ../src/flux_lib/mass/*.F90
$FC -c ../src/flux_lib/heat/*.F90
$FC -c ../src/flux_lib/momentum/*.F90
$FC -c ../src/flux_lib/*.F90
#$FC -shared -o flux_library.so ../src/flux_lib/*.F90
cd ..

87
src/.vscode/tasks.json vendored
Unescape View File

@ -0,0 +1,87 @@
{
// See https://go.microsoft.com/fwlink/?LinkId=733558
// for the documentation about the tasks.json format
"version": "2.0.0",
"tasks": [
{
"label": "build_flux_calculator",
"type": "shell",
"command": "cd /media/d/hagen/iow/climate/coupled_model/IOW_ESM/components/flux_calculator; rsync -r -i src mvkradtk@blogin:/scratch/usr/mvkradtk/IOW_ESM/components/flux_calculator/.; ssh -t mvkradtk@blogin 'cd /scratch/usr/mvkradtk/IOW_ESM/components/flux_calculator/; source build.sh'",
"group": {
"kind": "build",
"isDefault": false
},
"presentation": {
"reveal": "always",
"panel": "new"
},
//"problemMatcher": "$msCompile"
"problemMatcher": {
"owner": "ifort",
"fileLocation": ["relative", "${workspaceRoot}"],
"pattern": {
"regexp": "^../src/(.*?)\\((.*?)\\): (warning|error) (.*)",
"file": 1,
"line": 2,
"severity": 3,
"message": 4
}
}
},
{
"label": "build_flux_calculator_local",
"type": "shell",
"command": "cd /media/d/hagen/iow/climate/coupled_model/IOW_ESM/components/flux_calculator; source build_local.sh",
"group": {
"kind": "build",
"isDefault": false
},
"presentation": {
"reveal": "always",
"panel": "new"
},
"problemMatcher": "$msCompile"
},
{
"label": "build_flux_calculator_haumea",
"type": "shell",
"command": "cd /media/d/hagen/iow/climate/coupled_model/IOW_ESM/components/flux_calculator; rsync -r -i -u src haumea1:/data/hr275/IOW_ESM/components/flux_calculator/.; ssh -t haumea1 'cd /data/hr275/IOW_ESM/components/flux_calculator/; srun -t 01:00:00 -p compute -N1 --tasks-per-node 1 --pty bash -c \"source ~/.bash_profile; source build_debug.sh\"'",
"group": {
"kind": "build",
"isDefault": true
},
"presentation": {
"reveal": "always",
"panel": "new"
},
"problemMatcher": {
"owner": "ifort",
"fileLocation": [
"relative",
"${workspaceRoot}"
],
"pattern": {
"regexp": "^../src/(.*?)\\((.*?)\\): (warning|error) (.*)",
"file": 1,
"line": 2,
"severity": 3,
"message": 4
}
}
},
{
"label": "run_IOW_ESM_haumea",
"type": "shell",
"command": "ssh -t haumea1 'cd /data/hr275/IOW_ESM/scripts/run; source ~/.bash_profile; sbatch jobscript'",
"group": {
"kind": "test",
"isDefault": true
},
"presentation": {
"reveal": "always",
"panel": "new"
},
"problemMatcher": []
}
]
}

52
src/decomp_def.F90
Unescape View File

@ -0,0 +1,52 @@
!****************************************************************************************
SUBROUTINE decomp_def(id_paral,id_size,id_im,id_jm,id_rank,id_npes,id_unit)
!
IMPLICIT NONE
INTEGER, DIMENSION(id_size), INTENT(out) :: id_paral(id_size)
INTEGER, INTENT(in) :: id_size
INTEGER, INTENT(in) :: id_im ! Grid dimension in i
INTEGER, INTENT(in) :: id_jm ! Grid dimension in j
INTEGER, INTENT(in) :: id_rank ! Rank of process
INTEGER, INTENT(in) :: id_npes ! Number of processes involved in the coupling
INTEGER, INTENT(in) :: id_unit ! Unit of log file
INTEGER :: il_imjm, il_partj
!
il_imjm = id_im*id_jm
il_partj = id_jm/id_npes ! Nbr of latitude circles in the partition
!
#ifdef DECOMP_APPLE
! Each process is responsible for a part of field defined by
! the number of grid points and the offset of the first point
!
WRITE (id_unit,*) 'APPLE partitioning'
!
IF (id_rank .LT. (id_npes-1)) THEN
id_paral (1) = 1
id_paral (2) = id_rank*(il_partj * id_im)
id_paral (3) = il_partj * id_im
ELSE
id_paral (1) = 1
id_paral (2) = id_rank*(il_partj * id_im)
id_paral (3) = il_imjm-(id_rank*(il_partj * id_im))
ENDIF
!
#elif defined DECOMP_BOX
!
WRITE (id_unit,*) 'BOX partitioning'
!
! Each process is responsible for a rectangular box of il_partj lines
!
id_paral (1) = 2
id_paral (5) = id_im
id_paral (2) = id_rank*il_partj*id_im
id_paral (3) = id_im
IF (id_rank .LT. (id_npes-1)) THEN
id_paral (4) = il_partj
ELSE
id_paral (4) = id_jm-(id_rank*il_partj)
ENDIF
!
#endif
!
END SUBROUTINE decomp_def

1241
src/flux_calculator.F90
View File

File diff suppressed because it is too large Load Diff

560
src/flux_calculator_basic.F90
Unescape View File

@ -0,0 +1,560 @@
MODULE flux_calculator_basic
! This module contains:
! - type definitions
! - size limits for arrays
! - the list of variable names and their indexes
! - basic auxiliary functions
IMPLICIT NONE
!!!!!!!!!! FUNCTIONS DEFINED IN THIS MODULE
PUBLIC add_input_field
PUBLIC allocate_localvar
PUBLIC distribute_input_field
PUBLIC do_regridding
PUBLIC init_localvar
PUBLIC init_varname_idx
PUBLIC nullify_localvars
PUBLIC prepare_regridding
!!!!!!!!!! NOW EVERYTHING ELSE
#ifdef NO_USE_DOUBLE_PRECISION
INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(6,37) ! real
#elif defined USE_DOUBLE_PRECISION
INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307) ! double
#endif
INTEGER, PARAMETER :: MAX_BOTTOM_MODELS = 10 ! how many models we can couple to the atmospheric model
INTEGER, PARAMETER :: MAX_SURFACE_TYPES = 10 ! how many surface types (land or ice classes) are allowed
INTEGER, PARAMETER :: MAX_TASKS_PER_MODEL = 1000 ! on how many PEs per bottom model can the flux calculator run
INTEGER, PARAMETER :: MAX_VARS = 100 ! how many variables can a model send/receive on each grid
INTEGER, PARAMETER :: MAX_INPUT_FIELDS = 1000 ! how many variables can flux_calculator receive in total
INTEGER, PARAMETER :: MAX_OUTPUT_FIELDS = 1000 ! how many variables can flux_calculator send in total
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! If you add a varname ????, you need to do 4 changes: !
! (1) increase MAX_VARNAMES !
! (2) add '????' to varnames !
! (3) declare idx_varname !
! (4) add a line in init_varname_idx !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
INTEGER, PARAMETER :: MAX_VARNAMES = 30
CHARACTER(len=4), PARAMETER, DIMENSION(MAX_VARNAMES) :: varnames = [ &
'ALBE', 'AMOI', 'AMOM', 'FARE', 'FICE', 'PATM', 'PSUR', &
'QATM', 'TATM', 'TSUR', 'UATM', 'VATM', 'U10M', 'V10M', & ! variables read in
'QSUR', & ! auxiliary variables calculated
'HLAT', 'HSEN', & ! heat fluxes
'MEVA', 'MPRE', 'MRAI', 'MSNO', & ! mass fluxes
'RBBR', 'RLWD', 'RLWU', 'RSID', 'RSIU', 'RSIN', 'RSDD', & ! radiation fluxes, the last are: Shortwave_Indirect_Down, Shortwave_Indirect_Up, Shortwave_Indirect_Net, Shortwave_Direct_Down
'UMOM', 'VMOM'] ! momentum fluxes
INTEGER :: idx_ALBE, idx_AMOI, idx_AMOM, idx_FARE, idx_FICE, idx_PATM, idx_PSUR
INTEGER :: idx_QATM, idx_TATM, idx_TSUR, idx_UATM, idx_VATM, idx_U10M, idx_V10M
INTEGER :: idx_QSUR
INTEGER :: idx_HLAT, idx_HSEN
INTEGER :: idx_MEVA, idx_MPRE, idx_MRAI, idx_MSNO
INTEGER :: idx_RBBR, idx_RLWD, idx_RLWU, idx_RSID, idx_RSIU, idx_RSIN, idx_RSDD
INTEGER :: idx_UMOM, idx_VMOM
CHARACTER(len=2), DIMENSION(0:MAX_SURFACE_TYPES) :: numtype
CHARACTER(len=6), DIMENSION(3) :: grid_name = ['t_grid', 'u_grid', 'v_grid']
INTEGER :: w_unit ! a logfile to write the progress and error messages
TYPE integerarray
INTEGER (kind=4), DIMENSION(:), POINTER :: field
LOGICAL :: allocated
END TYPE integerarray
TYPE realarray
REAL (kind=wp), DIMENSION(:), POINTER :: field
LOGICAL :: allocated
LOGICAL :: put_to_t_grid ! whether this shall be regridded to the t_grid after reading/calculation here
LOGICAL :: put_to_u_grid ! whether this shall be regridded to the u_grid after reading/calculation here
LOGICAL :: put_to_v_grid ! whether this shall be regridded to the v_grid after reading/calculation here
END TYPE realarray
TYPE realarray2
REAL (kind=wp), DIMENSION(:,:), POINTER :: field
LOGICAL :: allocated
END TYPE realarray2
! define a type to store the variables
TYPE local_fields_type
! will be allocated with number of gridcells
TYPE(realarray), DIMENSION(MAX_VARNAMES) :: var ! these arrays will be allocated if the variable is present
END TYPE local_fields_type
! define a type to store the input/output fields
TYPE io_fields_type
REAL (kind=wp), DIMENSION(:), POINTER :: field ! arrays will just be pointers to local_fields arrays
CHARACTER(len=8) :: name ! will store the OASIS variable name
INTEGER :: which_grid ! 1=t_grid, 2=u_grid, 3=v_grid
INTEGER :: id ! the id returned by oasis_def_var
LOGICAL :: early ! whether or not we communicate these variable before all others
INTEGER :: surface_type ! surface type, used for regridding of input data
INTEGER :: idx ! idx_???? number of variable, used for regridding of input data
END TYPE io_fields_type
! define a type to store a regridding matrix (e.g. t_grid to u_grid)
TYPE sparse_regridding_matrix
INTEGER :: num_elements
TYPE(integerarray) :: src_index
TYPE(integerarray) :: dst_index
TYPE(realarray) :: weight
END TYPE sparse_regridding_matrix
CONTAINS
SUBROUTINE add_input_field(myname, my_letter, surface_type, which_grid, local_field, num_input_fields, input_field)
CHARACTER(len=4), INTENT(IN) :: myname ! e.g. "PATM"
CHARACTER(len=1), INTENT(IN) :: my_letter ! 'A' for atmosphere, my_bottom_letter for bottom model
INTEGER, INTENT(IN) :: surface_type
INTEGER, INTENT(IN) :: which_grid ! 1=t_grid, 2=u_grid, 3=v_grid
TYPE(local_fields_type), INTENT(INOUT) :: local_field ! here we want to set the pointer to
INTEGER, INTENT(INOUT) :: num_input_fields
type(io_fields_type), DIMENSION(:), INTENT(INOUT) :: input_field
INTEGER :: i
LOGICAL :: found
CHARACTER(len=2) :: appendstring ! two-digit number of surface type
appendstring = numtype(surface_type)
found=.FALSE.
DO i=1,MAX_VARNAMES
IF (varnames(i) == myname) THEN
found=.TRUE.
num_input_fields = num_input_fields + 1
IF (num_input_fields > MAX_INPUT_FIELDS) THEN
WRITE (w_unit,*) "Maximum number of input fields exceeded. Consider increasing MAX_INPUT_FIELDS and then recompile flux_calculator."
CALL mpi_finalize(1)
ELSE
input_field(num_input_fields)%name="R"//my_letter//myname//appendstring
input_field(num_input_fields)%which_grid = which_grid
input_field(num_input_fields)%field => local_field%var(i)%field
input_field(num_input_fields)%early=.FALSE.
IF ((myname=="FARE") .OR. (myname=="TSUR")) THEN
input_field(num_input_fields)%early=.TRUE.
ENDIF
input_field(num_input_fields)%surface_type = surface_type
input_field(num_input_fields)%idx = i
ENDIF
ENDIF
ENDDO
IF (.NOT. found) THEN
WRITE (w_unit,*) "Could not add input field for variable ",myname," because flux_calculator does not know this variable."
CALL mpi_finalize(1)
ENDIF
END SUBROUTINE add_input_field
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE add_output_field(myname, my_letter, surface_type, which_grid, num_surface_types, grid_size, &
uniform, default_value, local_field, num_input_fields, input_field, num_output_fields, output_field)
CHARACTER(len=4), INTENT(IN) :: myname ! e.g. "PATM"
CHARACTER(len=1), INTENT(IN) :: my_letter ! "A" for atmospheric model, my_bottom_letter for bottom model
INTEGER, INTENT(IN) :: surface_type ! 0=entire cell, 1,2,3,... = specific surface type
INTEGER, INTENT(IN) :: which_grid ! 1=t_grid, 2=u_grid, 3=v_grid
INTEGER, INTENT(IN) :: num_surface_types ! how many does my bottom model have
INTEGER, INTENT(IN) :: grid_size
LOGICAL, INTENT(IN) :: uniform ! whether this flux is the same for each surface_type
REAL(kind=wp), INTENT(IN) :: default_value ! if no method calculated the flux
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field ! here we want to set the pointer to
INTEGER, INTENT(IN) :: num_input_fields
type(io_fields_type), DIMENSION(:), INTENT(IN) :: input_field
INTEGER, INTENT(INOUT) :: num_output_fields
type(io_fields_type), DIMENSION(:), INTENT(INOUT) :: output_field
INTEGER :: i,j
LOGICAL :: found, found_all_fluxes, found_all_areas
CHARACTER(len=2) :: appendstring ! two-digit number of surface type
appendstring=numtype(surface_type)
! First check if this already exists as input field. In this case we do not need the output field.
found=.FALSE.
DO i=1,num_input_fields
IF (trim(input_field(i)%name)=="R"//my_letter//myname//appendstring) found=.TRUE.
ENDDO
IF (.NOT. found) THEN ! Okay field does not come as input
DO i=1,MAX_VARNAMES
IF (varnames(i) == myname) THEN
found=.TRUE.
! Do we want output for surface_type=0 (the entire exchange grid cell)?
IF (surface_type==0) THEN
! Yes, output for the entire cell
IF (uniform) THEN
! Flux is the same everywhere so we only need to find one surface_type for which it is calculated
DO j=1,num_surface_types
IF (ASSOCIATED(local_field(j,which_grid)%var(i)%field) .AND. (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field))) THEN
local_field(0,which_grid)%var(i)%field => local_field(j,which_grid)%var(i)%field
ENDIF
ENDDO
ELSE
! Flux will differ, so we will need a value for each surface_type and an area fraction of that surface_type
found_all_fluxes = .TRUE.
found_all_areas = .TRUE.
DO j=1,num_surface_types
IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(i)%field)) found_all_fluxes=.FALSE.
IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(idx_FARE)%field)) found_all_areas=.FALSE.
ENDDO
IF (.NOT. found_all_areas) THEN
! fractional area is not provided - this is fatal
WRITE (w_unit,*) "ERROR: Output field ",myname," has not been defined as uniform (flux_?_uniform=.FALSE.)."
WRITE (w_unit,*) "To calculate its average value across different surface_types, their fractional area (FARE) must be given but is missing."
CALL mpi_finalize(1)
ELSE
IF (found_all_fluxes) THEN
! this is nice, we found everything we need, so just allocate
IF (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field)) THEN
ALLOCATE(local_field(0,which_grid)%var(i)%field(grid_size))
local_field(0,which_grid)%var(i)%allocated=.TRUE.
ENDIF
ENDIF
ENDIF
ENDIF
IF (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field)) THEN
! Okay we have not found a pointer yet or allocated the field we want to send.
! But we also have not encountered a serious error.
! So we may just be missing the calculation of this flux.
! This will just give a warning message because the fluxes have standard values (typically zero)
WRITE (w_unit,*) "WARNING: Flux ",myname," cannot be calculated for surface_type=0 => set to ",default_value
ALLOCATE(local_field(0,which_grid)%var(i)%field(grid_size))
local_field(0,which_grid)%var(i)%field = default_value
local_field(0,which_grid)%var(i)%allocated=.TRUE.
ENDIF
ELSE
IF (uniform) THEN
IF (.NOT. ASSOCIATED(local_field(surface_type,which_grid)%var(i)%field)) THEN
!okay it is not calculated for this surface_type, but probably for another surface_type, so we can use its values
DO j=1,num_surface_types
IF (ASSOCIATED(local_field(j,which_grid)%var(i)%field) .AND. (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field))) THEN
local_field(0,which_grid)%var(i)%field => local_field(j,which_grid)%var(i)%field
ENDIF
ENDDO
ENDIF
ENDIF
! We want output for one specific surface type only
IF (.NOT. ASSOCIATED(local_field(surface_type,which_grid)%var(i)%field)) THEN
! We are missing the calculation of this flux.
! This will just give a warning message because the fluxes have standard values (typically zero)
WRITE (w_unit,*) "WARNING: Flux ",myname," cannot be calculated for surface_type=",surface_type," => set to ",default_value
ALLOCATE(local_field(surface_type,which_grid)%var(i)%field(grid_size))
local_field(surface_type,which_grid)%var(i)%field = default_value
local_field(surface_type,which_grid)%var(i)%allocated=.TRUE.
ENDIF
ENDIF
num_output_fields = num_output_fields + 1
IF (num_output_fields > MAX_OUTPUT_FIELDS) THEN
WRITE (w_unit,*) "Maximum number of output fields exceeded. Consider increasing MAX_OUTPUT_FIELDS and then recompile flux_calculator."
CALL mpi_finalize(1)
ELSE
output_field(num_output_fields)%name="S"//my_letter//myname//appendstring
output_field(num_output_fields)%which_grid = which_grid
output_field(num_output_fields)%field => local_field(surface_type,which_grid)%var(i)%field
output_field(num_output_fields)%early=.FALSE.
IF ((myname=="RBBR") .OR. (myname=="TSUR")) THEN
output_field(num_output_fields)%early=.TRUE.
ENDIF
ENDIF
ENDIF
ENDDO
IF (.NOT. found) THEN
WRITE (w_unit,*) "Could not add output field for variable ",myname," because flux_calculator does not know this variable."
CALL mpi_finalize(1)
ENDIF
ENDIF
END SUBROUTINE add_output_field
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE allocate_localvar(myname, length, local_field)
CHARACTER(len=4), INTENT(IN) :: myname ! e.g. "PATM"
INTEGER, INTENT(IN) :: length
TYPE(local_fields_type), INTENT(INOUT) :: local_field
INTEGER :: i
LOGICAL :: found
found=.FALSE.
DO i=1,MAX_VARNAMES
IF (varnames(i) == myname) THEN
found=.TRUE.
ALLOCATE(local_field%var(i)%field(length))
local_field%var(i)%allocated = .TRUE.
local_field%var(i)%put_to_t_grid = .FALSE.
local_field%var(i)%put_to_u_grid = .FALSE.
local_field%var(i)%put_to_v_grid = .FALSE.
ENDIF
ENDDO
IF (.NOT. found) THEN
WRITE (w_unit,*) "Could not allocate local variable ",myname," because flux_calculator does not know this variable."
CALL mpi_finalize(1)
ENDIF
END SUBROUTINE allocate_localvar
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_localvar(myname, value, local_field)
CHARACTER(len=4), INTENT(IN) :: myname ! e.g. "PATM"
REAL(kind=wp), INTENT(IN) :: value
TYPE(local_fields_type), INTENT(INOUT) :: local_field
INTEGER :: i
LOGICAL :: found
found=.FALSE.
DO i=1,MAX_VARNAMES
IF (varnames(i) == myname) THEN
found=.TRUE.
local_field%var(i)%field = value
ENDIF
ENDDO
IF (.NOT. found) THEN
WRITE (w_unit,*) "Could not set local variable ",myname," to constant value because flux_calculator does not know this variable."
CALL mpi_finalize(1)
ENDIF
END SUBROUTINE init_localvar
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE distribute_input_field(myname, which_grid, from_surface_type, to_surface_type, num_surface_types, local_field)
CHARACTER(len=4), INTENT(IN) :: myname ! e.g. "PATM"
INTEGER, INTENT(IN) :: which_grid
INTEGER, INTENT(IN) :: from_surface_type
INTEGER, INTENT(IN) :: to_surface_type
INTEGER, INTENT(IN) :: num_surface_types
TYPE(local_fields_type), INTENT(INOUT), DIMENSION(0:,:) :: local_field
INTEGER :: i, j
LOGICAL :: found
found=.FALSE.
DO i=1,MAX_VARNAMES
IF (varnames(i) == myname) THEN
found=.TRUE.
DO j=1,num_surface_types
IF ((j == to_surface_type) .OR. ((j /= from_surface_type) .AND. (to_surface_type == 0))) THEN
local_field(j,which_grid)%var(i)%field => local_field(from_surface_type,which_grid)%var(i)%field
ENDIF
ENDDO
ENDIF
ENDDO
IF (.NOT. found) THEN
WRITE (w_unit,*) "Could not distribute local variable ",myname," to other surface_types because flux_calculator does not know this variable."
CALL mpi_finalize(1)
ENDIF
END SUBROUTINE distribute_input_field
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE prepare_regridding(varidx, surface_type, local_field, my_bottom_model, &
regrid_u_to_t, regrid_v_to_t, regrid_t_to_u, regrid_t_to_v, &
grid_size)
INTEGER, INTENT(IN) :: varidx
INTEGER, INTENT(IN) :: surface_type ! set to 0 for all surface types
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
INTEGER, INTENT(IN) :: my_bottom_model
CHARACTER(len=4), DIMENSION(:,:,:), INTENT(IN) :: regrid_u_to_t
CHARACTER(len=4), DIMENSION(:,:,:), INTENT(IN) :: regrid_v_to_t
CHARACTER(len=4), DIMENSION(:,:,:), INTENT(IN) :: regrid_t_to_u
CHARACTER(len=4), DIMENSION(:,:,:), INTENT(IN) :: regrid_t_to_v
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
INTEGER :: i, j, k
INTEGER :: from_grid, to_grid
i=varidx
from_grid=2; to_grid=1
DO j=1,MAX_SURFACE_TYPES
IF ((j==surface_type) .OR. (surface_type==0)) THEN
DO k=1,MAX_VARS
IF (varnames(i) == trim(regrid_u_to_t(my_bottom_model,j,k))) THEN
IF (local_field(j,to_grid)%var(i)%allocated) THEN
! This should not be the case - the variable already exists on the destination grid
WRITE (w_unit,*) "Could not regrid local variable ",varnames(i)," from ",grid_name(from_grid)," to ",grid_name(to_grid),&
" as requested in the namelist, because it already exists on that grid."
CALL mpi_finalize(1)
ELSE
ALLOCATE(local_field(j,to_grid)%var(i)%field(grid_size(to_grid)))
local_field(j,to_grid)%var(i)%allocated = .TRUE.
local_field(j,from_grid)%var(i)%put_to_t_grid = .TRUE.
WRITE(w_unit,*) ' u_grid -> t_grid: ',varnames(i),' for surface type ',j
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
from_grid=3; to_grid=1
DO j=1,MAX_SURFACE_TYPES
IF ((j==surface_type) .OR. (surface_type==0)) THEN
DO k=1,MAX_VARS
IF (varnames(i) == trim(regrid_v_to_t(my_bottom_model,j,k))) THEN
IF (local_field(j,to_grid)%var(i)%allocated) THEN
! This should not be the case - the variable already exists on the destination grid
WRITE (w_unit,*) "Could not regrid local variable ",varnames(i)," from ",grid_name(from_grid)," to ",grid_name(to_grid),&
" as requested in the namelist, because it already exists on that grid."
CALL mpi_finalize(1)
ELSE
ALLOCATE(local_field(j,to_grid)%var(i)%field(grid_size(to_grid)))
local_field(j,to_grid)%var(i)%allocated = .TRUE.
local_field(j,from_grid)%var(i)%put_to_t_grid = .TRUE.
WRITE(w_unit,*) ' v_grid -> t_grid: ',varnames(i),' for surface type ',j
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
from_grid=1; to_grid=2
DO j=1,MAX_SURFACE_TYPES
IF ((j==surface_type) .OR. (surface_type==0)) THEN
DO k=1,MAX_VARS
IF (varnames(i) == trim(regrid_t_to_u(my_bottom_model,j,k))) THEN
IF (local_field(j,to_grid)%var(i)%allocated) THEN
! This should not be the case - the variable already exists on the destination grid
WRITE (w_unit,*) "Could not regrid local variable ",varnames(i)," from ",grid_name(from_grid)," to ",grid_name(to_grid),&
" as requested in the namelist, because it already exists on that grid."
CALL mpi_finalize(1)
ELSE
ALLOCATE(local_field(j,to_grid)%var(i)%field(grid_size(to_grid)))
local_field(j,to_grid)%var(i)%allocated = .TRUE.
local_field(j,from_grid)%var(i)%put_to_u_grid = .TRUE.
WRITE(w_unit,*) ' t_grid -> u_grid: ',varnames(i),' for surface type ',j
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
from_grid=1; to_grid=3
DO j=1,MAX_SURFACE_TYPES
IF ((j==surface_type) .OR. (surface_type==0)) THEN
DO k=1,MAX_VARS
IF (varnames(i) == trim(regrid_t_to_v(my_bottom_model,j,k))) THEN
IF (local_field(j,to_grid)%var(i)%allocated) THEN
! This should not be the case - the variable already exists on the destination grid
WRITE (w_unit,*) "Could not regrid local variable ",varnames(i)," from ",grid_name(from_grid)," to ",grid_name(to_grid),&
" as requested in the namelist, because it already exists on that grid."
CALL mpi_finalize(1)
ELSE
ALLOCATE(local_field(j,to_grid)%var(i)%field(grid_size(to_grid)))
local_field(j,to_grid)%var(i)%allocated = .TRUE.
local_field(j,from_grid)%var(i)%put_to_v_grid = .TRUE.
WRITE(w_unit,*) ' t_grid -> v_grid: ',varnames(i),' for surface type ',j
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
END SUBROUTINE prepare_regridding
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE do_regridding(varidx, surface_type, local_field, &
regrid_u_to_t_matrix, regrid_v_to_t_matrix, regrid_t_to_u_matrix, regrid_t_to_v_matrix)
INTEGER, INTENT(IN) :: varidx
INTEGER, INTENT(IN) :: surface_type ! 0 stands for all surface types
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
TYPE(sparse_regridding_matrix), INTENT(IN) :: regrid_u_to_t_matrix
TYPE(sparse_regridding_matrix), INTENT(IN) :: regrid_v_to_t_matrix
TYPE(sparse_regridding_matrix), INTENT(IN) :: regrid_t_to_u_matrix
TYPE(sparse_regridding_matrix), INTENT(IN) :: regrid_t_to_v_matrix
INTEGER :: i, j, k, from, to
INTEGER :: from_grid, to_grid
i=varidx
DO j=1,MAX_SURFACE_TYPES
IF ((j==surface_type) .OR. (surface_type==0)) THEN
IF (local_field(j,2)%var(i)%put_to_t_grid) THEN
from=2; to=1
local_field(j,to)%var(i)%field = 0.0
DO k=1,regrid_u_to_t_matrix%num_elements
local_field(j,to)%var(i)%field(regrid_u_to_t_matrix%dst_index%field(k)) = & ! sparse matrix multiplication
local_field(j,to)%var(i)%field(regrid_u_to_t_matrix%dst_index%field(k)) + &
local_field(j,from)%var(i)%field(regrid_u_to_t_matrix%src_index%field(k)) * &
regrid_u_to_t_matrix%weight%field(k)
ENDDO
ENDIF
IF (local_field(j,3)%var(i)%put_to_t_grid) THEN
from=3; to=1
local_field(j,to)%var(i)%field = 0.0
DO k=1,regrid_v_to_t_matrix%num_elements
local_field(j,to)%var(i)%field(regrid_v_to_t_matrix%dst_index%field(k)) = & ! sparse matrix multiplication
local_field(j,to)%var(i)%field(regrid_v_to_t_matrix%dst_index%field(k)) + &
local_field(j,from)%var(i)%field(regrid_v_to_t_matrix%src_index%field(k)) * &
regrid_v_to_t_matrix%weight%field(k)
ENDDO
ENDIF
IF (local_field(j,1)%var(i)%put_to_u_grid) THEN
from=1; to=2
local_field(j,to)%var(i)%field = 0.0
DO k=1,regrid_t_to_u_matrix%num_elements
local_field(j,to)%var(i)%field(regrid_t_to_u_matrix%dst_index%field(k)) = & ! sparse matrix multiplication
local_field(j,to)%var(i)%field(regrid_t_to_u_matrix%dst_index%field(k)) + &
local_field(j,from)%var(i)%field(regrid_t_to_u_matrix%src_index%field(k)) * &
regrid_t_to_u_matrix%weight%field(k)
ENDDO
ENDIF
IF (local_field(j,1)%var(i)%put_to_v_grid) THEN
from=1; to=3
local_field(j,to)%var(i)%field = 0.0
DO k=1,regrid_t_to_v_matrix%num_elements
local_field(j,to)%var(i)%field(regrid_t_to_v_matrix%dst_index%field(k)) = & ! sparse matrix multiplication
local_field(j,to)%var(i)%field(regrid_t_to_v_matrix%dst_index%field(k)) + &
local_field(j,from)%var(i)%field(regrid_t_to_v_matrix%src_index%field(k)) * &
regrid_t_to_v_matrix%weight%field(k)
ENDDO
ENDIF
ENDIF
ENDDO
END SUBROUTINE do_regridding
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE init_varname_idx
INTEGER :: i
DO i=1,MAX_VARNAMES
! INPUT VARIABLES:
IF (varnames(i) == 'ALBE') idx_ALBE = i ! Albedo of the surface (1)
IF (varnames(i) == 'AMOI') idx_AMOI = i ! Diffusion coefficient for moisture (1)
IF (varnames(i) == 'AMOM') idx_AMOM = i ! Diffusion coefficient for momentum (1)
IF (varnames(i) == 'FARE') idx_FARE = i ! Fraction of grid cell area covered by this surface_type (1)
IF (varnames(i) == 'FICE') idx_FICE = i ! Fraction of area covered by ice (1)
IF (varnames(i) == 'PATM') idx_PATM = i ! Pressure in lowest atmospheric grid cell (Pa)
IF (varnames(i) == 'PSUR') idx_PSUR = i ! Pressure at surface (Pa)
IF (varnames(i) == 'QATM') idx_QATM = i ! Moisture content of lowest atmospheric grid cell (kg/kg)
IF (varnames(i) == 'TATM') idx_TATM = i ! Temperature in lowest atmospheric grid cell (K)
IF (varnames(i) == 'TSUR') idx_TSUR = i ! Temperature at the surface (K)
IF (varnames(i) == 'UATM') idx_UATM = i ! Eastward velocity in lowest atmospheric grid cell (m/s)
IF (varnames(i) == 'VATM') idx_VATM = i ! Northward velocity in lowest atmospheric grid cell (m/s)
IF (varnames(i) == 'U10M') idx_U10M = i ! Eastward velocity 10m above the surface (m/s)
IF (varnames(i) == 'V10M') idx_V10M = i ! Northward velocity 10m above the surface (m/s)
! AUXILIARY VARIABLES:
IF (varnames(i) == 'QSUR') idx_QSUR = i ! Moisture content directly above the surface (kg/kg)
! FLUXES: (all are positive upward, such that e.g. precipitation is always negative)
IF (varnames(i) == 'HLAT') idx_HLAT = i ! Heat flux (latent) (W/m2)
IF (varnames(i) == 'HSEN') idx_HSEN = i ! Heat flux (sensible) (W/m2)
IF (varnames(i) == 'MEVA') idx_MEVA = i ! Mass flux (evaporation minus condensation) (kg/m2/s)
IF (varnames(i) == 'MPRE') idx_MPRE = i ! Mass flux (total precipitation, negative values) (kg/m2/s)
IF (varnames(i) == 'MRAI') idx_MRAI = i ! Mass flux (liquid precipitation "rain", neg. val.) (kg/m2/s)
IF (varnames(i) == 'MSNO') idx_MSNO = i ! Mass flux (solid precipitation "snow", neg. val.) (kg/m2/s)
IF (varnames(i) == 'RBBR') idx_RBBR = i ! Radiation flux (blackbody radiation of surface) (W/m2)
IF (varnames(i) == 'RLWD') idx_RLWD = i ! Radiation flux (longwave downward, negative values) (W/m2)
IF (varnames(i) == 'RLWU') idx_RLWU = i ! Radiation flux (longwave upward) (W/m2)
IF (varnames(i) == 'RSID') idx_RSID = i ! Radiation flux (shortwave indirect downward, neg. val.) (W/m2)
IF (varnames(i) == 'RSIU') idx_RSIU = i ! Radiation flux (shortwave indirect upward) (W/m2)
IF (varnames(i) == 'RSIN') idx_RSIN = i ! Radiation flux (shortwave indirect net upward) (W/m2)
IF (varnames(i) == 'RSDD') idx_RSID = i ! Radiation flux (shortwave directed downward, neg. val.) (W/m2)
IF (varnames(i) == 'UMOM') idx_UMOM = i ! Upward flux of eastward momentum (N/m2)
IF (varnames(i) == 'VMOM') idx_VMOM = i ! Upward flux of northward momentum (N/m2)
ENDDO
END SUBROUTINE init_varname_idx
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE nullify_localvars(local_field)
TYPE(local_fields_type), INTENT(INOUT) :: local_field
INTEGER :: i
DO i=1,MAX_VARNAMES
NULLIFY(local_field%var(i)%field) ! this ensures that ASSOCIATED(local_field%var(i)%field)==.FALSE.
local_field%var(i)%allocated = .FALSE.
local_field%var(i)%put_to_t_grid = .FALSE.
local_field%var(i)%put_to_u_grid = .FALSE.
local_field%var(i)%put_to_v_grid = .FALSE.
ENDDO
END SUBROUTINE nullify_localvars
END MODULE flux_calculator_basic

272
src/flux_calculator_calculate.F90
Unescape View File

@ -0,0 +1,272 @@
MODULE flux_calculator_calculate
! This module contains functions that do the calculations using the flux library functions
use flux_calculator_basic
use flux_library
IMPLICIT NONE
!!!!!!!!!! FUNCTIONS DEFINED IN THIS MODULE
PUBLIC calc_spec_vapor_surface
PUBLIC calc_flux_mass_evap
PUBLIC calc_flux_radiation_blackbody
PUBLIC average_across_surface_types
!!!!!!!!!! NOW EVERYTHING ELSE
CONTAINS
!!!!!!!!!! AUXILIARY VARIABLES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE calc_spec_vapor_surface(my_bottom_model, num_surface_types, which_grid, methods, grid_size, local_field)
! calculates specific water vapor content on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
INTEGER, INTENT(IN) :: which_grid
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'QSUR'
CHARACTER(len=20) :: method
INTEGER :: i, j
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='CCLM') THEN
DO j=1,grid_size(which_grid)
CALL spec_vapor_surface_cclm(local_field(i,which_grid)%var(idx_QSUR)%field(j), &
local_field(i,which_grid)%var(idx_FICE)%field(j), &
local_field(i,which_grid)%var(idx_PSUR)%field(j), &
local_field(i,which_grid)%var(idx_TSUR)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
END SUBROUTINE calc_spec_vapor_surface
!!!!!!!!!! MASS FLUXES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE calc_flux_mass_evap(my_bottom_model, num_surface_types, methods, grid_size, local_field)
! calculates specific water vapor content on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'MEVA'
CHARACTER(len=20) :: method
INTEGER :: i, j
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,1)%var(idx_MEVA)%field(:) = 0.0
ELSEIF (trim(method)=='CCLM') THEN
DO j=1,grid_size(1)
CALL flux_mass_evap_cclm(local_field(i,1)%var(idx_MEVA)%field(j), &
local_field(i,1)%var(idx_AMOI)%field(j), &
local_field(i,1)%var(idx_PSUR)%field(j), &
local_field(i,1)%var(idx_QATM)%field(j), &
local_field(i,1)%var(idx_QSUR)%field(j), &
local_field(i,1)%var(idx_TATM)%field(j), &
local_field(i,1)%var(idx_UATM)%field(j), &
local_field(i,1)%var(idx_VATM)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_MEVA,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_mass_evap
!!!!!!!!!! HEAT FLUXES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE calc_flux_heat_latent(my_bottom_model, num_surface_types, methods, grid_size, local_field)
! calculates latent heat flux on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'HLAT'
CHARACTER(len=20) :: method
INTEGER :: i, j
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,1)%var(idx_MEVA)%field(:) = 0.0
ELSEIF (trim(method)=='water') THEN
DO j=1,grid_size(1)
CALL flux_heat_latent_water(local_field(i,1)%var(idx_HLAT)%field(j), &
local_field(i,1)%var(idx_MEVA)%field(j))
ENDDO
ELSEIF (trim(method)=='ice') THEN
DO j=1,grid_size(1)
CALL flux_heat_latent_ice(local_field(i,1)%var(idx_HLAT)%field(j), &
local_field(i,1)%var(idx_MEVA)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_HLAT,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_heat_latent
SUBROUTINE calc_flux_heat_sensible(my_bottom_model, num_surface_types, methods, grid_size, local_field)
! calculates sensible heat flux on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'HSEN'
CHARACTER(len=20) :: method
INTEGER :: i, j
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,1)%var(idx_MEVA)%field(:) = 0.0
ELSEIF (trim(method)=='CCLM') THEN
DO j=1,grid_size(1)
CALL flux_heat_sensible_cclm(local_field(i,1)%var(idx_HSEN)%field(j), &
local_field(i,1)%var(idx_AMOI)%field(j), &
local_field(i,1)%var(idx_PATM)%field(j), &
local_field(i,1)%var(idx_PSUR)%field(j), &
local_field(i,1)%var(idx_QATM)%field(j), &
local_field(i,1)%var(idx_TATM)%field(j), &
local_field(i,1)%var(idx_TSUR)%field(j), &
local_field(i,1)%var(idx_UATM)%field(j), &
local_field(i,1)%var(idx_VATM)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_HSEN,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_heat_sensible
!!!!!!!!!! MOMENTUM FLUXES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE calc_flux_momentum_east(my_bottom_model, num_surface_types, which_grid, methods, grid_size, local_field)
! calculates sensible heat flux on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
INTEGER, INTENT(IN) :: which_grid
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'UMOM'
CHARACTER(len=20) :: method
INTEGER :: i, j
REAL :: dummy ! since northward momentum may be required on another grid
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,which_grid)%var(idx_UMOM)%field(:) = 0.0
ELSEIF (trim(method)=='CCLM') THEN
DO j=1,grid_size(which_grid)
CALL flux_momentum_cclm(local_field(i,which_grid)%var(idx_UMOM)%field(j), &
dummy, &
local_field(i,which_grid)%var(idx_AMOM)%field(j), &
local_field(i,which_grid)%var(idx_PSUR)%field(j), &
local_field(i,which_grid)%var(idx_QSUR)%field(j), &
local_field(i,which_grid)%var(idx_TSUR)%field(j), &
local_field(i,which_grid)%var(idx_UATM)%field(j), &
local_field(i,which_grid)%var(idx_VATM)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_UMOM,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_momentum_east
SUBROUTINE calc_flux_momentum_north(my_bottom_model, num_surface_types, which_grid, methods, grid_size, local_field)
! calculates sensible heat flux on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
INTEGER, INTENT(IN) :: which_grid
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'VMOM'
CHARACTER(len=20) :: method
INTEGER :: i, j
REAL :: dummy ! since northward momentum may be required on another grid
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,which_grid)%var(idx_VMOM)%field(:) = 0.0
ELSEIF (trim(method)=='CCLM') THEN
DO j=1,grid_size(which_grid)
CALL flux_momentum_cclm(dummy, &
local_field(i,which_grid)%var(idx_VMOM)%field(j), &
local_field(i,which_grid)%var(idx_AMOM)%field(j), &
local_field(i,which_grid)%var(idx_PSUR)%field(j), &
local_field(i,which_grid)%var(idx_QSUR)%field(j), &
local_field(i,which_grid)%var(idx_TSUR)%field(j), &
local_field(i,which_grid)%var(idx_UATM)%field(j), &
local_field(i,which_grid)%var(idx_VATM)%field(j))
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_VMOM,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_momentum_north
!!!!!!!!!! RADIATION FLUXES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE calc_flux_radiation_blackbody(my_bottom_model, num_surface_types, methods, grid_size, local_field)
! calculates blackbody radiation on t_grid
INTEGER, INTENT(IN) :: my_bottom_model
INTEGER, INTENT(IN) :: num_surface_types
CHARACTER(len=20), DIMENSION(:,:), INTENT(IN) :: methods
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
CHARACTER(len=4), PARAMETER :: myvarname = 'RBBR'
CHARACTER(len=20) :: method
INTEGER :: i, j
DO i=1,num_surface_types
method = methods(my_bottom_model,i)
IF (trim(method) /= 'none') THEN
IF (trim(method)=='zero') THEN
local_field(i,1)%var(idx_RBBR)%field(:) = 0.0
ELSEIF (trim(method)=='StBo') THEN
DO j=1,grid_size(1)
CALL flux_radiation_blackbody_StBo(local_field(i,1)%var(idx_RBBR)%field(j), &
local_field(i,1)%var(idx_TSUR)%field(j) )
ENDDO
ENDIF
ENDIF
ENDDO
CALL average_across_surface_types(1,idx_RBBR,num_surface_types,grid_size,local_field)
END SUBROUTINE calc_flux_radiation_blackbody
!!!!!!!!!! AVERAGING ROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE average_across_surface_types(which_grid, my_idx, num_surface_types, grid_size, local_field)
INTEGER, INTENT(IN) :: which_grid
INTEGER, INTENT(IN) :: my_idx
INTEGER, INTENT(IN) :: num_surface_types
INTEGER, DIMENSION(:), INTENT(IN) :: grid_size
TYPE(local_fields_type), DIMENSION(0:,:), INTENT(INOUT) :: local_field
INTEGER :: i,j
IF (local_field(0,1)%var(my_idx)%allocated) THEN
local_field(0,1)%var(my_idx)%field=0.0
DO i=1,num_surface_types
DO j=1,grid_size(which_grid)
local_field(0,1)%var(my_idx)%field(j) = local_field(0,1)%var(my_idx)%field(j) + &
local_field(i,1)%var(my_idx)%field(j)*local_field(i,1)%var(idx_FARE)%field(j)
ENDDO
ENDDO
ENDIF
END SUBROUTINE average_across_surface_types
END MODULE flux_calculator_calculate

169
src/flux_calculator_io.F90
Unescape View File

@ -0,0 +1,169 @@
MODULE flux_calculator_io
! This module contains functions for NetCDF IO which flux_calculator performs
use flux_calculator_basic
IMPLICIT NONE
public read_scrip_grid_dimensions, read_regridding_matrix
contains
SUBROUTINE read_scrip_grid_dimensions(grid_filename, mype, npes, num_bottom_models, num_tasks_per_model, &
num_grid_cells, grid_size_global, grid_size, grid_offset, grid_lon, grid_lat, grid_area)
! Read a SCRIP-type NetCDF exchange grid file and return info on which gridcells we need to calculate
! in this MPI instance.
USE netcdf
CHARACTER(len=50), INTENT(IN) :: grid_filename
INTEGER, INTENT(IN) :: mype
INTEGER, INTENT(IN) :: npes
INTEGER, INTENT(IN) :: num_bottom_models
INTEGER, DIMENSION(:), INTENT(IN) :: num_tasks_per_model
INTEGER, DIMENSION(:,:), INTENT(INOUT) :: num_grid_cells ! this works in two ways:
! (a) num_bottom_models=1 and num_tasks_per_model=1, then this is an output field and its value (all points of this grid) is assigned here
! (b) otherwise this is an input field (that needs to be defined in the namelist and will be passed here)
INTEGER, INTENT(OUT) :: grid_size_global ! grid size for the entire exchange grid (all bottom models)
INTEGER, INTENT(OUT) :: grid_size ! grid size for this task only
INTEGER, INTENT(OUT) :: grid_offset ! sum of grid sizes from all preceding tasks
! (each grid cell of the exchange grid is used by a single task only, so when one bottom model's exchange grid is split to different tasks,
! some cells may need to be duplicated in the exchange grid file due to necessary regridding (u,v <-> t) from the "halos")
TYPE(realarray2), INTENT(INOUT) :: grid_lon
TYPE(realarray2), INTENT(INOUT) :: grid_lat
TYPE(realarray2), INTENT(INOUT) :: grid_area
INTEGER :: nc ! NetCDF file id
INTEGER :: dimid_grid_size ! NetCDF dimension id
INTEGER :: varid ! NetCDF variable id
CALL hdlerr(NF90_OPEN(grid_filename, NF90_NOWRITE, nc), __LINE__ )
CALL hdlerr( NF90_INQ_DIMID(nc, 'grid_size' , dimid_grid_size), __LINE__ ) ! get variable id
CALL hdlerr( NF90_INQUIRE_DIMENSION(ncid=nc,dimid=dimid_grid_size,len=grid_size_global), __LINE__ )
ALLOCATE(grid_lon%field(grid_size_global,1)); grid_lon%allocated = .TRUE.
ALLOCATE(grid_lat%field(grid_size_global,1)); grid_lat%allocated = .TRUE.
ALLOCATE(grid_area%field(grid_size_global,1)); grid_area%allocated = .TRUE.
CALL hdlerr( NF90_INQ_VARID(nc, 'grid_center_lon' , varid), __LINE__ ) ! get variable id
CALL hdlerr( NF90_GET_VAR (nc, varid, grid_lon%field(:,1), [1], [grid_size_global]), __LINE__ ) ! get variable values
CALL hdlerr( NF90_INQ_VARID(nc, 'grid_center_lat' , varid), __LINE__ )
CALL hdlerr( NF90_GET_VAR (nc, varid, grid_lat%field(:,1), [1], [grid_size_global]), __LINE__ )
CALL hdlerr( NF90_INQ_VARID(nc, 'grid_area' , varid), __LINE__ )
CALL hdlerr( NF90_GET_VAR (nc, varid, grid_area%field(:,1), [1], [grid_size_global]), __LINE__ )
CALL hdlerr(NF90_CLOSE(nc), __LINE__ )
! most simple case: single model, single task -> use full grid
IF (num_bottom_models == 1) THEN
IF (num_tasks_per_model(1) == 1) THEN
num_grid_cells = 0
num_grid_cells(1,1) = grid_size_global
grid_size = grid_size_global
grid_offset = 0
ENDIF
ENDIF
END SUBROUTINE read_scrip_grid_dimensions
SUBROUTINE read_regridding_matrix(regridding_filename, &
src_grid_size, src_grid_offset, &
dst_grid_size, dst_grid_offset, &
mymatrix)
USE netcdf
INCLUDE 'mpif.h'
! read the sparse regridding matrix (in coordinates format) for the full exchange grid and extract
! those elements which are for the current task
CHARACTER(len=50), INTENT(IN) :: regridding_filename
INTEGER, INTENT(IN) :: src_grid_size
INTEGER, INTENT(IN) :: src_grid_offset
INTEGER, INTENT(IN) :: dst_grid_size
INTEGER, INTENT(IN) :: dst_grid_offset
TYPE(sparse_regridding_matrix), INTENT(INOUT) :: mymatrix
INTEGER :: i, num_total, num_matching
INTEGER :: nc ! NetCDF file id
INTEGER :: dimid_num_links ! NetCDF dimension id
INTEGER :: varid ! NetCDF variable id
TYPE(integerarray) :: all_src_index
TYPE(integerarray) :: all_dst_index
TYPE(realarray2) :: all_weight
! read in the overall exchange grid remapping matrix from the file
CALL hdlerr(NF90_OPEN(regridding_filename, NF90_NOWRITE, nc), __LINE__ )
CALL hdlerr( NF90_INQ_DIMID(nc, 'num_links' , dimid_num_links), __LINE__ ) ! get dimension id
CALL hdlerr( NF90_INQUIRE_DIMENSION(ncid=nc,dimid=dimid_num_links,len=num_total), __LINE__ )
ALLOCATE(all_src_index%field(num_total)); all_src_index%allocated = .TRUE.
ALLOCATE(all_dst_index%field(num_total)); all_dst_index%allocated = .TRUE.
ALLOCATE(all_weight%field(1,num_total)); all_weight%allocated = .TRUE.
CALL hdlerr( NF90_INQ_VARID(nc, 'src_address' , varid), __LINE__ ) ! get variable id
CALL hdlerr( NF90_GET_VAR (nc, varid, all_src_index%field(:), [1], [num_total]), __LINE__ ) ! get variable values
CALL hdlerr( NF90_INQ_VARID(nc, 'dst_address' , varid), __LINE__ ) ! get variable id
CALL hdlerr( NF90_GET_VAR (nc, varid, all_dst_index%field(:), [1], [num_total]), __LINE__ ) ! get variable values
CALL hdlerr( NF90_INQ_VARID(nc, 'remap_matrix' , varid), __LINE__ ) ! get variable id
CALL hdlerr( NF90_GET_VAR (nc, varid, all_weight%field(:,:), [1,1], [1,num_total]), __LINE__ ) ! get variable values
CALL hdlerr(NF90_CLOSE(nc), __LINE__ )
! find out how many entries are matching, so we can allocate the final array to the correct size
num_matching = 0
DO i=1,num_total
IF ((all_dst_index%field(i) .gt. dst_grid_offset) .AND. &
(all_dst_index%field(i) .le. dst_grid_offset + dst_grid_size)) THEN
num_matching = num_matching + 1
END IF
END DO
mymatrix%num_elements = num_matching
ALLOCATE(mymatrix%src_index%field(num_matching)); mymatrix%src_index%allocated = .TRUE.
ALLOCATE(mymatrix%dst_index%field(num_matching)); mymatrix%dst_index%allocated = .TRUE.
ALLOCATE(mymatrix%weight%field(num_matching)); mymatrix%weight%allocated = .TRUE.
! walk through the full matrix and copy only the matching entries to the matrix that will be used in this task
num_matching = 0
DO i=1,num_total
IF ((all_dst_index%field(i) .gt. dst_grid_offset) .AND. &
(all_dst_index%field(i) .le. dst_grid_offset + dst_grid_size)) THEN
num_matching = num_matching + 1
mymatrix%dst_index%field(num_matching) = all_dst_index%field(i)
mymatrix%src_index%field(num_matching) = all_src_index%field(i)
mymatrix%weight%field(num_matching) = all_weight%field(1,i)
END IF
END DO
! correct for the offset and check if all source cells are present at this PE
DO i=1,num_matching
mymatrix%dst_index%field(i) = mymatrix%dst_index%field(i) - dst_grid_offset
mymatrix%src_index%field(i) = mymatrix%src_index%field(i) - src_grid_offset
IF ((mymatrix%src_index%field(i) .LT. 1) .OR. (mymatrix%src_index%field(i) .GT. src_grid_size)) THEN
write ( * , * ) 'Regridding matrix did not match task decomposition in file ',regridding_filename
write ( * , * ) 'Stopped '
call MPI_Abort ( MPI_COMM_WORLD, 1, 1 )
END IF
END DO
! clean up: deallocate full arrays
DEALLOCATE(all_src_index%field)
DEALLOCATE(all_dst_index%field)
DEALLOCATE(all_weight%field)
END