Skip to content

Commit 4867f73

Browse files
fix: vtec extrapolation calculation
- result of 4 corners was not being stored correctly, overwritting corner11 - domain of ZG was incorrect - units of ZG for need to be m in the vtec calculation fixes #851
1 parent 4eea462 commit 4867f73

File tree

1 file changed

+5
-5
lines changed

1 file changed

+5
-5
lines changed

models/tiegcm/model_mod.f90

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77
! - Nick Dietrich fix_mmr. When do to this?
88
! - model_time
99
! - get_state_meta_data 2D variables
10-
! - test vtec
1110

1211
module model_mod
1312

@@ -379,9 +378,9 @@ subroutine model_interpolate(state_handle, ens_size, location, iqty, obs_val, is
379378
if ( iqty == QTY_VERTICAL_TEC ) then ! extrapolate vtec
380379

381380
call extrapolate_vtec(state_handle, ens_size, lon_below, lat_below, val11)
382-
call extrapolate_vtec(state_handle, ens_size, lon_below, lat_above, val11)
383-
call extrapolate_vtec(state_handle, ens_size, lon_above, lat_below, val11)
384-
call extrapolate_vtec(state_handle, ens_size, lon_above, lat_above, val11)
381+
call extrapolate_vtec(state_handle, ens_size, lon_below, lat_above, val12)
382+
call extrapolate_vtec(state_handle, ens_size, lon_above, lat_below, val21)
383+
call extrapolate_vtec(state_handle, ens_size, lon_above, lat_above, val22)
385384
obs_val(:) = interpolate(ens_size, lon_fract, lat_fract, val11, val12, val21, val22)
386385
istatus(:) = 0
387386

@@ -1111,7 +1110,7 @@ subroutine extrapolate_vtec(state_handle, ens_size, lon_index, lat_index, vTEC)
11111110
! ZG (interfaces)
11121111
do i = 1, nilev
11131112
idx = get_dart_vector_index(lon_index,lat_index, i, &
1114-
domain_id(RESTART_DOM), var_id)
1113+
domain_id(SECONDARY_DOM), ivarZG)
11151114
ZG(i, :) = get_state(idx, state_handle)
11161115
enddo
11171116

@@ -1135,6 +1134,7 @@ subroutine extrapolate_vtec(state_handle, ens_size, lon_index, lat_index, vTEC)
11351134

11361135
earth_radiusm = earth_radius * 1000.0_r8 ! Convert earth_radius in km to m
11371136
NE = NE * 1.0e+6_r8 ! Convert NE in #/cm^3 to #/m^3
1137+
ZG = ZG * 1.0e-2_r8 ! Convert ZG in cm to m
11381138

11391139
! Gravity at the top layer
11401140
GRAVITYtop(:) = gravity * (earth_radiusm / (earth_radiusm + ZG(nilev,:))) ** 2

0 commit comments

Comments
 (0)