FV3-JEDI
surface_variables_mod.f90
Go to the documentation of this file.
1 ! (C) Copyright 2018-2019 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 module surface_vt_mod
7 
8 use fv3jedi_geom_mod, only: fv3jedi_geom
9 use fv3jedi_kinds_mod, only: kind_real
10 use crtm_module, only: crtm_irlandcoeff_classification
11 use fv3jedi_constants_mod, only: rad2deg, deg2rad, pi
12 
13 use unstructured_interpolation_mod, only: unstrc_interp
14 
15 implicit none
16 private
17 
18 public crtm_surface
19 
20 contains
21 
22 !----------------------------------------------------------------------------
23 ! Surface quantities in the form needed by the crtm -------------------------
24 !----------------------------------------------------------------------------
25 
26 subroutine crtm_surface( geom, field_slmsk, field_sheleg, field_tsea, field_vtype, field_stype, &
27  field_vfrac, field_stc, field_smc, field_u_srf, field_v_srf, field_f10m, &
28  field_sss, &
29  land_type, vegetation_type, soil_type, water_coverage, land_coverage, &
30  ice_coverage, snow_coverage, lai, water_temperature, land_temperature, &
31  ice_temperature, snow_temperature, soil_moisture_content, &
32  vegetation_fraction, soil_temperature, snow_depth, wind_speed, &
33  wind_direction, sea_surface_salinity)
34 
35 implicit none
36 
37 !Arguments
38 type(fv3jedi_geom) , intent(in) :: geom
39 real(kind=kind_real), intent(in) :: field_slmsk(geom%isc:geom%iec,geom%jsc:geom%jec,1)
40 real(kind=kind_real), intent(in) :: field_sheleg(geom%isc:geom%iec,geom%jsc:geom%jec,1)
41 real(kind=kind_real), intent(in) :: field_tsea(geom%isc:geom%iec,geom%jsc:geom%jec,1)
42 real(kind=kind_real), intent(in) :: field_vtype(geom%isc:geom%iec,geom%jsc:geom%jec,1)
43 real(kind=kind_real), intent(in) :: field_stype(geom%isc:geom%iec,geom%jsc:geom%jec,1)
44 real(kind=kind_real), intent(in) :: field_vfrac(geom%isc:geom%iec,geom%jsc:geom%jec,1)
45 real(kind=kind_real), intent(in) :: field_stc(geom%isc:geom%iec,geom%jsc:geom%jec,1)
46 real(kind=kind_real), intent(in) :: field_smc(geom%isc:geom%iec,geom%jsc:geom%jec,1)
47 real(kind=kind_real), intent(in) :: field_u_srf(geom%isc:geom%iec,geom%jsc:geom%jec,1)
48 real(kind=kind_real), intent(in) :: field_v_srf(geom%isc:geom%iec,geom%jsc:geom%jec,1)
49 real(kind=kind_real), intent(in) :: field_f10m(geom%isc:geom%iec,geom%jsc:geom%jec,1)
50 real(kind=kind_real), intent(in) :: field_sss(geom%isc:geom%iec,geom%jsc:geom%jec,1)
51 real(kind=kind_real), intent(out) :: vegetation_type(geom%isc:geom%iec,geom%jsc:geom%jec,1)
52 real(kind=kind_real), intent(out) :: land_type(geom%isc:geom%iec,geom%jsc:geom%jec,1)
53 real(kind=kind_real), intent(out) :: soil_type(geom%isc:geom%iec,geom%jsc:geom%jec,1)
54 real(kind=kind_real), intent(out) :: water_coverage(geom%isc:geom%iec,geom%jsc:geom%jec,1)
55 real(kind=kind_real), intent(out) :: land_coverage(geom%isc:geom%iec,geom%jsc:geom%jec,1)
56 real(kind=kind_real), intent(out) :: ice_coverage(geom%isc:geom%iec,geom%jsc:geom%jec,1)
57 real(kind=kind_real), intent(out) :: snow_coverage(geom%isc:geom%iec,geom%jsc:geom%jec,1)
58 real(kind=kind_real), intent(out) :: lai(geom%isc:geom%iec,geom%jsc:geom%jec,1)
59 real(kind=kind_real), intent(out) :: water_temperature(geom%isc:geom%iec,geom%jsc:geom%jec,1)
60 real(kind=kind_real), intent(out) :: land_temperature(geom%isc:geom%iec,geom%jsc:geom%jec,1)
61 real(kind=kind_real), intent(out) :: ice_temperature(geom%isc:geom%iec,geom%jsc:geom%jec,1)
62 real(kind=kind_real), intent(out) :: snow_temperature(geom%isc:geom%iec,geom%jsc:geom%jec,1)
63 real(kind=kind_real), intent(out) :: soil_moisture_content(geom%isc:geom%iec,geom%jsc:geom%jec,1)
64 real(kind=kind_real), intent(out) :: vegetation_fraction(geom%isc:geom%iec,geom%jsc:geom%jec,1)
65 real(kind=kind_real), intent(out) :: soil_temperature(geom%isc:geom%iec,geom%jsc:geom%jec,1)
66 real(kind=kind_real), intent(out) :: snow_depth(geom%isc:geom%iec,geom%jsc:geom%jec,1)
67 real(kind=kind_real), intent(out) :: wind_speed(geom%isc:geom%iec,geom%jsc:geom%jec,1)
68 real(kind=kind_real), intent(out) :: wind_direction(geom%isc:geom%iec,geom%jsc:geom%jec,1)
69 real(kind=kind_real), intent(out) :: sea_surface_salinity(geom%isc:geom%iec,geom%jsc:geom%jec,1)
70 
71 !Locals
72 real(kind=kind_real), parameter :: minsnow = 1.0_kind_real / 10.0_kind_real
73 real(kind=kind_real), parameter :: windlimit = 0.0001_kind_real
74 real(kind=kind_real), parameter :: quadcof(4, 2 ) = &
75  reshape((/0.0_kind_real, 1.0_kind_real, 1.0_kind_real, 2.0_kind_real, &
76  1.0_kind_real, -1.0_kind_real, 1.0_kind_real, -1.0_kind_real/), (/4, 2/))
77 
78 integer :: itype, istype
79 integer :: istyp00
80 integer :: isflg
81 integer :: lai_type, iquadrant
82 logical :: lwind
83 real(kind=kind_real) :: sfcpct(0:3), ts(0:3), wgtavg(0:3), dtskin(0:3)
84 real(kind=kind_real) :: sno00
85 real(kind=kind_real) :: sst00
86 real(kind=kind_real) :: ss00
87 real(kind=kind_real) :: tsavg,ssavg
88 real(kind=kind_real) :: vty, sty, vfr, stp, sm, sn, ss
89 real(kind=kind_real) :: uu5, vv5, f10, sfc_speed, windratio, windangle, windscale
90 real(kind=kind_real) :: wind10, wind10_direction
91 
92 !From GSI
93 integer, parameter :: gfs_soil_n_types = 9
94 integer, parameter :: gfs_vegetation_n_types = 13
95 integer, parameter :: invalid_land = 0
96 integer, parameter :: compacted_soil = 1
97 integer, parameter :: tilled_soil = 2
98 integer, parameter :: irrigated_low_vegetation = 5
99 integer, parameter :: meadow_grass = 6
100 integer, parameter :: scrub = 7
101 integer, parameter :: broadleaf_forest = 8
102 integer, parameter :: pine_forest = 9
103 integer, parameter :: tundra = 10
104 integer, parameter :: grass_soil = 11
105 integer, parameter :: broadleaf_pine_forest = 12
106 integer, parameter :: grass_scrub = 13
107 integer, parameter :: urban_concrete = 15
108 integer, parameter :: broadleaf_brush = 17
109 integer, parameter :: wet_soil = 18
110 integer, parameter :: scrub_soil = 19
111 integer, parameter :: nvege_type = 20
112 integer, parameter :: igbp_n_types = 20
113 integer, parameter :: soil_n_types = 16
114 integer, allocatable,dimension(:) :: map_to_crtm_ir
115 integer, allocatable,dimension(:) :: map_to_crtm_mwave
116 integer, parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_gfs=(/4, &
117  1, 5, 2, 3, 8, 9, 6, 6, 7, 8, 12, 7, 12, 13, 11, 0, 10, 10, 11/)
118  integer, parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_npoess=(/pine_forest, &
119  broadleaf_forest, pine_forest, broadleaf_forest, broadleaf_pine_forest, &
120  scrub, scrub_soil, broadleaf_brush, broadleaf_brush, scrub, broadleaf_brush, &
121  tilled_soil, urban_concrete, tilled_soil, invalid_land, compacted_soil, &
122  invalid_land, tundra, tundra, tundra/)
123  integer, parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_igbp=(/1, &
124  2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, &
125  20/)
126  integer, parameter, dimension(1:SOIL_N_TYPES) :: map_soil_to_crtm=(/1, &
127  1, 4, 2, 2, 8, 7, 2, 6, 5, 2, 3, 8, 1, 6, 9/)
128 
129 integer :: ji, jj
130 integer :: slmsk
131 integer :: vtype
132 integer :: stype
133 real(kind=kind_real) :: sheleg
134 real(kind=kind_real) :: tsea
135 real(kind=kind_real) :: vfrac
136 real(kind=kind_real) :: snwdph
137 real(kind=kind_real) :: stc
138 real(kind=kind_real) :: smc
139 real(kind=kind_real) :: u_srf
140 real(kind=kind_real) :: v_srf
141 real(kind=kind_real) :: f10m
142 real(kind=kind_real) :: sss
143 
144 ! Zero outputs
145 vegetation_type = 0.0_kind_real
146 land_type = 0.0_kind_real
147 soil_type = 0.0_kind_real
148 water_coverage = 0.0_kind_real
149 land_coverage = 0.0_kind_real
150 ice_coverage = 0.0_kind_real
151 snow_coverage = 0.0_kind_real
152 lai = 0.0_kind_real
153 water_temperature = 0.0_kind_real
154 land_temperature = 0.0_kind_real
155 ice_temperature = 0.0_kind_real
156 snow_temperature = 0.0_kind_real
157 soil_moisture_content = 0.0_kind_real
158 vegetation_fraction = 0.0_kind_real
159 soil_temperature = 0.0_kind_real
160 snow_depth = 0.0_kind_real
161 wind_speed = 0.0_kind_real
162 wind_direction = 0.0_kind_real
163 sea_surface_salinity = 0.0_kind_real
164 
165 ! Vegetation maps
166 allocate(map_to_crtm_ir(nvege_type))
167 allocate(map_to_crtm_mwave(nvege_type))
168 map_to_crtm_ir = igbp_to_igbp
169 map_to_crtm_mwave = igbp_to_gfs
170 !TODO, this belongs in ufo or with advanced locations
171 !select case ( TRIM(CRTM_IRlandCoeff_Classification()) )
172 ! case('NPOESS'); map_to_crtm_ir=igbp_to_npoess
173 ! case('IGBP') ; map_to_crtm_ir=igbp_to_igbp
174 !end select
175 
176 ! Loop over all grid points
177 do jj = geom%jsc, geom%jec
178  do ji = geom%isc, geom%iec
179 
180  slmsk = nint(field_slmsk(ji,jj,1))
181  vtype = nint(field_vtype(ji,jj,1))
182  stype = nint(field_stype(ji,jj,1))
183  sheleg = field_sheleg(ji,jj,1)
184  tsea = field_tsea(ji,jj,1)
185  vfrac = field_vfrac(ji,jj,1)
186  stc = field_stc(ji,jj,1)
187  smc = field_smc(ji,jj,1)
188  u_srf = field_u_srf(ji,jj,1)
189  v_srf = field_v_srf(ji,jj,1)
190  f10m = field_f10m(ji,jj,1)
191  sss = field_sss(ji,jj,1)
192 
193  dtskin = 0.0_kind_real !TODO need real skin temperature increment?
194 
195  lwind = .true.
196 
197  ! Stage 1, like deter_sfc in GSI
198  ! ------------------------------
199  istyp00 = slmsk
200  sno00 = sheleg !sno00 = snwdph
201  sst00 = tsea
202  tsavg = sst00
203  ss00 = sss
204 
205  ssavg = ss00
206 
207  if (istyp00 >=1 .and. sno00 > minsnow) istyp00 = 3
208 
209  sfcpct = 0.0_kind_real
210  sfcpct(istyp00) = 1.0
211 
212  isflg = 0
213  if(sfcpct(0) > 0.99_kind_real)then
214  isflg = 0
215  else if(sfcpct(1) > 0.99_kind_real)then
216  isflg = 1
217  else if(sfcpct(2) > 0.99_kind_real)then
218  isflg = 2
219  else if(sfcpct(3) > 0.99_kind_real)then
220  isflg = 3
221  else
222  isflg = 4
223  end if
224 
225  ts(0:3)=0.0_kind_real
226  wgtavg(0:3)=0.0_kind_real
227  vfr=0.0_kind_real
228  stp=0.0_kind_real
229  sty=0.0_kind_real
230  vty=0.0_kind_real
231  sm=0.0_kind_real
232  sn=0.0_kind_real
233  ss=0.0_kind_real
234 
235  if(istyp00 == 1)then
236  vty = vtype
237  sty = stype
238  wgtavg(1) = wgtavg(1) + 1.0
239  ts(1)=ts(1)+sst00
240  vfr =vfr + vfrac
241  stp =stp + stc
242  sm =sm + smc
243  else if(istyp00 == 2)then
244  wgtavg(2) = wgtavg(2) + 1.0
245  ts(2)=ts(2)+sst00
246  else if(istyp00 == 3)then
247  wgtavg(3) = wgtavg(3) + 1.0
248  ts(3)=ts(3)+sst00
249  sn = sn + sno00
250  else
251  wgtavg(0) = wgtavg(0) + 1.0
252  ts(0)=ts(0)+sst00
253  ss =ss +ss00
254  end if
255 
256  if(wgtavg(0) > 0.0_kind_real)then
257  ts(0) = ts(0)/wgtavg(0)
258  ss = ss/wgtavg(0)
259  else
260  ts(0) = tsavg
261  ss = ssavg
262  end if
263 
264  if(wgtavg(1) > 0.0_kind_real)then
265  ts(1) = ts(1)/wgtavg(1)
266  sm = sm/wgtavg(1)
267  vfr = vfr/wgtavg(1)
268  stp = stp/wgtavg(1)
269  else
270  ts(1) = tsavg
271  sm=1.0_kind_real
272  end if
273 
274  if(wgtavg(2) > 0.0_kind_real)then
275  ts(2) = ts(2)/wgtavg(2)
276  else
277  ts(2) = tsavg
278  end if
279 
280  if(wgtavg(3) > 0.0_kind_real)then
281  ts(3) = ts(3)/wgtavg(3)
282  sn = sn/wgtavg(3)
283  else
284  ts(3) = tsavg
285  end if
286 
287  f10 = f10m
288 
289  ! Stage 2 - like crtm_interface from GSI
290  ! --------------------------------------
291  ! If vty/sty will give maximum values, that will be out of range for CRTM, set to 1
292  if (vty == 15) vty = 1
293  if (sty == 16) sty = 1
294 
295  itype = vty
296  istype = sty
297 
298  itype = min(max(1,itype),nvege_type)
299  istype = min(max(1,istype),soil_n_types)
300  land_type(ji,jj,1) = real(max(1,map_to_crtm_mwave(itype)),kind_real)
301  vegetation_type(ji,jj,1) = real(max(1,map_to_crtm_mwave(itype)),kind_real)
302  soil_type(ji,jj,1) = real(map_soil_to_crtm(istype),kind_real)
303  lai_type = real(map_to_crtm_mwave(itype),kind_real)
304 
305  water_coverage(ji,jj,1) = min(max(0.0_kind_real,sfcpct(0)),1.0_kind_real)
306  land_coverage(ji,jj,1) = min(max(0.0_kind_real,sfcpct(1)),1.0_kind_real)
307  ice_coverage(ji,jj,1) = min(max(0.0_kind_real,sfcpct(2)),1.0_kind_real)
308  snow_coverage(ji,jj,1) = min(max(0.0_kind_real,sfcpct(3)),1.0_kind_real)
309 
310  lai(ji,jj,1) = 0.0_kind_real
311 
312  if (land_coverage(ji,jj,1) > 0.0_kind_real) then
313 
314  if(lai_type>0)then
315  call get_lai(lai_type,lai(ji,jj,1)) !TODO: does nothing yet
316  endif
317 
318  ! for Glacial land ice soil type and vegetation type
319  if(soil_type(ji,jj,1) == 9 .OR. vegetation_type(ji,jj,1) == 13) then
320  ice_coverage(ji,jj,1) = min(ice_coverage(ji,jj,1) + land_coverage(ji,jj,1), 1.0_kind_real)
321  land_coverage(ji,jj,1) = 0.0_kind_real
322  endif
323 
324  endif
325 
326  if (lwind) then
327 
328  !Interpolate lowest level winds to observation location
329  uu5 = u_srf
330  vv5 = v_srf
331 
332  sfc_speed = f10*sqrt(uu5*uu5+vv5*vv5)
333  wind10 = sfc_speed
334  if (uu5*f10 >= 0.0_kind_real .and. vv5*f10 >= 0.0_kind_real) iquadrant = 1
335  if (uu5*f10 >= 0.0_kind_real .and. vv5*f10 < 0.0_kind_real) iquadrant = 2
336  if (uu5*f10 < 0.0_kind_real .and. vv5*f10 >= 0.0_kind_real) iquadrant = 4
337  if (uu5*f10 < 0.0_kind_real .and. vv5*f10 < 0.0_kind_real) iquadrant = 3
338  if (abs(vv5*f10) >= windlimit) then
339  windratio = (uu5*f10) / (vv5*f10)
340  else
341  windratio = 0.0_kind_real
342  if (abs(uu5*f10) > windlimit) then
343  windratio = windscale * uu5*f10
344  endif
345  endif
346  windangle = atan(abs(windratio)) ! wind azimuth is in radians
347  wind10_direction = quadcof(iquadrant, 1) * pi + windangle * quadcof(iquadrant, 2)
348  wind_speed(ji,jj,1) = sfc_speed
349  wind_direction(ji,jj,1) = rad2deg*wind10_direction
350 
351  else
352 
353  wind_speed(ji,jj,1) = 0.0_kind_real
354  wind_direction(ji,jj,1) = 0.0_kind_real
355 
356  endif
357 
358  water_temperature(ji,jj,1) = max(ts(0) + dtskin(0), 270._kind_real)
359  sea_surface_salinity(ji,jj,1) = ss
360 
361  !TODO, is nst_gsi ever > 1?
362  !if(nst_gsi > 1 .and. water_coverage(1) > 0.0_kind_real) then
363  !water_temperature(ji,jj,1) = max(data_s(itref)+data_s(idtw)-data_s(idtc) + dtskin(0), 271._kind_real)
364  !endif
365 
366  land_temperature(ji,jj,1) = ts(1) + dtskin(1)
367  ice_temperature(ji,jj,1) = min(ts(2) + dtskin(2), 280._kind_real)
368  snow_temperature(ji,jj,1) = min(ts(3) + dtskin(3), 280._kind_real)
369  soil_moisture_content(ji,jj,1) = sm
370  vegetation_fraction(ji,jj,1) = vfr
371  soil_temperature(ji,jj,1) = stp
372  snow_depth(ji,jj,1) = sn
373 
374  enddo
375 enddo
376 
377 deallocate(map_to_crtm_ir)
378 deallocate(map_to_crtm_mwave)
379 
380 end subroutine crtm_surface
381 
382 !----------------------------------------------------------------------------
383 
384 subroutine get_lai(lai_type,lai)
385 
386  implicit none
387 
388  integer , intent(in ) :: lai_type
389  real(kind=kind_real), intent(out) :: lai
390 
391  !Dummy code, needs to be figured out
392  if (lai_type .ne. 0) then
393  lai = 0.0_kind_real
394  endif
395 
396  end subroutine get_lai
397 
398 !----------------------------------------------------------------------------
399 
400 end module surface_vt_mod
surface_vt_mod
Definition: surface_variables_mod.f90:6
fv3jedi_constants_mod::rad2deg
real(kind=kind_real), parameter, public rad2deg
Definition: fv3jedi_constants_mod.f90:13
fv3jedi_geom_mod
Fortran module handling geometry for the FV3 model.
Definition: fv3jedi_geom_mod.f90:8
fv3jedi_geom_mod::fv3jedi_geom
Fortran derived type to hold geometry data for the FV3JEDI model.
Definition: fv3jedi_geom_mod.f90:46
fv3jedi_constants_mod
Definition: fv3jedi_constants_mod.f90:6
surface_vt_mod::crtm_surface
subroutine, public crtm_surface(geom, field_slmsk, field_sheleg, field_tsea, field_vtype, field_stype, field_vfrac, field_stc, field_smc, field_u_srf, field_v_srf, field_f10m, field_sss, land_type, vegetation_type, soil_type, water_coverage, land_coverage, ice_coverage, snow_coverage, lai, water_temperature, land_temperature, ice_temperature, snow_temperature, soil_moisture_content, vegetation_fraction, soil_temperature, snow_depth, wind_speed, wind_direction, sea_surface_salinity)
Definition: surface_variables_mod.f90:34
fv3jedi_constants_mod::deg2rad
real(kind=kind_real), parameter, public deg2rad
Definition: fv3jedi_constants_mod.f90:14
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:14
fv3jedi_kinds_mod
Definition: fv3jedi_kinds_mod.f90:6
fv3jedi_constants_mod::pi
real(kind=kind_real), parameter, public pi
Definition: fv3jedi_constants_mod.f90:16
surface_vt_mod::get_lai
subroutine get_lai(lai_type, lai)
Definition: surface_variables_mod.f90:385