1module oceanObservations_mod
2 ! MODULE oceanObservations_mod (prefix='oobs' category='1. High-level functionality')
3 !
4 !:Purpose: storage for ocean observations related subroutines.
5 !
6 use midasMpi_mod
7 use utilities_mod
8 use obsSpaceData_mod
9 use codtyp_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use horizontalCoord_mod
13 use verticalCoord_mod
14 use oceanMask_mod
15 use timeCoord_mod
16 use codePrecision_mod
17 use sqliteRead_mod
18 use bufr_mod
19
20 implicit none
21
22 save
23 private
24
25 ! Public functions/subroutines
26 public :: oobs_pseudoSST
27
28 ! External functions
29 integer, external :: fnom, fclos
30
31 ! mpi topology
32 integer :: myLatBeg, myLatEnd
33 integer :: myLonBeg, myLonEnd
34 integer :: latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
35
36 contains
37
38 !----------------------------------------------------------------------------------------
39 ! oobs_pseudoSST
40 !----------------------------------------------------------------------------------------
41 subroutine oobs_pseudoSST(hco, vco, iceFractionThreshold, outputSST, outputFreshWaterST, &
42 iceThinning, outputFileName, seaWaterThreshold)
43 !
44 !:Purpose: to generate pseudo SST data
45 !
46 implicit none
47
48 ! Arguments:
49 type(struct_hco), pointer , intent(inout) :: hco ! horizontal grid structure
50 type(struct_vco), pointer , intent(in) :: vco ! vertical grid structure
51 real(4) , intent(in) :: iceFractionThreshold ! consider no ice condition below this threshold
52 real(4) , intent(in) :: outputSST ! output SST value for pseudo observations
53 real(4) , intent(in) :: outputFreshWaterST ! output fresh water surface temperature for pseudo obs
54 integer , intent(in) :: iceThinning ! generate pseudo obs in every 'iceThinning' points
55 character(len=*) , intent(in) :: outputFileName
56 real(4) , intent(in) :: seaWaterThreshold ! to distinguish inland water from sea water
57
58 ! Locals:
59 type(struct_gsv) :: stateVector_ice, stateVector_seaWater
60 real(4), pointer :: seaIce_ptr(:, :, :), seaWater_ptr(:, :, :)
61 type(struct_ocm) :: oceanMask
62 integer :: numberIceCoveredPoints, lonIndex, latIndex, dateStamp, inlandWaterPoints
63 integer :: datePrint, timePrint, imode, seaWaterPoints
64 integer :: randomSeed, newDate, ierr
65 integer, allocatable :: iceDomainIndexesAux(:), iceDomainIndexes(:)
66 real(4), allocatable :: seaWaterFractionAux(:), iceLonsAux(:), iceLatsAux(:)
67 real(4), allocatable :: seaWaterFraction(:), iceLons(:), iceLats(:)
68 type(struct_obs) :: obsData
69
70 ! get mpi topology
71 call mmpi_setup_lonbands(hco%ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
72 call mmpi_setup_latbands(hco%nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
73
74 ! get latest sea-ice analysis
75 call gsv_allocate(stateVector_ice, 1, hco, vco, dataKind_opt = 4, &
76 datestamp_opt = -1, mpi_local_opt = .false., varNames_opt = (/'LG'/))
77 call gio_readFromFile(stateVector_ice, './seaice_analysis', ' ','A', &
78 unitConversion_opt=.false., containsFullField_opt=.true.)
79 call gsv_getField(stateVector_ice, seaIce_ptr)
80
81 ! read sea water fraction
82 call gsv_allocate(stateVector_seaWater, 1, hco, vco, dataKind_opt = 4, &
83 datestamp_opt = -1, mpi_local_opt = .false., varNames_opt = (/'VF'/))
84 call gio_readFromFile(stateVector_seaWater, './seaice_analysis', ' ','A', &
85 unitConversion_opt=.false., containsFullField_opt=.true.)
86 call gsv_getField(stateVector_seaWater, seaWater_ptr)
87
88 ! Get land mask from analysisgrid file (1=water, 0=land)
89 call ocm_readMaskFromFile(oceanMask, hco, vco, './analysisgrid')
90
91 allocate(iceDomainIndexesAux((myLonEnd - myLonBeg + 1) * (myLatEnd - myLatBeg + 1)))
92 allocate(iceLonsAux((myLonEnd - myLonBeg + 1) * (myLatEnd - myLatBeg + 1)))
93 allocate(iceLatsAux((myLonEnd - myLonBeg + 1) * (myLatEnd - myLatBeg + 1)))
94 allocate(seaWaterFractionAux((myLonEnd - myLonBeg + 1) * (myLatEnd - myLatBeg + 1)))
95
96 numberIceCoveredPoints = 0
97 inlandWaterPoints = 0
98 seaWaterPoints = 0
99
100 do lonIndex = myLonBeg, myLonEnd
101 do latIndex = myLatBeg, myLatEnd
102 if (oceanMask%mask(lonIndex, latIndex, 1)) then
103 if (seaice_ptr(lonIndex, latIndex, 1) > iceFractionThreshold) then
104
105 numberIceCoveredPoints = numberIceCoveredPoints + 1
106 iceDomainIndexesAux(numberIceCoveredPoints) = numberIceCoveredPoints
107 seaWaterFractionAux(numberIceCoveredPoints) = seaWater_ptr(lonIndex, latIndex, 1)
108
109 if(seaWater_ptr(lonIndex, latIndex, 1) <= seaWaterThreshold) then
110 inlandWaterPoints = inlandWaterPoints + 1
111 else
112 seaWaterPoints = seaWaterPoints + 1
113 end if
114
115 iceLonsAux(numberIceCoveredPoints) = hco%lon2d_4 (lonIndex, latIndex)
116 iceLatsAux(numberIceCoveredPoints) = hco%lat2d_4 (lonIndex, latIndex)
117
118 end if
119 end if
120 end do
121 end do
122 call ocm_deallocate(oceanMask)
123 call gsv_deallocate(stateVector_ice)
124 call gsv_deallocate(stateVector_seaWater)
125 write(*,*) 'oobs_pseudoSST: ', numberIceCoveredPoints, ' ice-covered points found'
126 write(*,*) 'oobs_pseudoSST: where ', inlandWaterPoints, ' are inland water points'
127 write(*,*) 'oobs_pseudoSST: ', seaWaterPoints, ' sea water points found'
128
129 if (numberIceCoveredPoints > 0) then
130 allocate(iceDomainIndexes(1:numberIceCoveredPoints))
131 iceDomainIndexes(:) = iceDomainIndexesAux(1:numberIceCoveredPoints)
132 allocate(seaWaterFraction(1:numberIceCoveredPoints))
133 seaWaterFraction(:) = seaWaterFractionAux(1:numberIceCoveredPoints)
134 allocate(iceLons(1:numberIceCoveredPoints))
135 allocate(iceLats(1:numberIceCoveredPoints))
136 iceLons(:) = iceLonsAux(1:numberIceCoveredPoints)
137 iceLats(:) = iceLatsAux(1:numberIceCoveredPoints)
138 end if
139 deallocate(iceLonsAux)
140 deallocate(iceLatsAux)
141 deallocate(iceDomainIndexesAux)
142 deallocate(seaWaterFractionAux)
143
144 dateStamp = tim_getDatestampFromFile('./seaice_analysis', varNameForDate_opt = 'LG')
145 write(*,*) 'oobs_pseudoSST: datestamp: ', dateStamp
146 ! compute random seed from the date for randomly forming sea-ice subdomain
147 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
148 ierr = newdate(dateStamp, datePrint, timePrint, imode)
149 timePrint = timePrint / 1000000
150 datePrint = datePrint * 100 + timePrint
151
152 ! Remove the century, keeping 2 digits of the year
153 randomSeed = datePrint - 100000000 * (datePrint / 100000000)
154 write(*,*) 'oobs_pseudoSST: datePrint, timePrint: ', datePrint, timePrint
155
156 if (numberIceCoveredPoints > 0) then
157
158 call utl_randomOrderInt(iceDomainIndexes, randomSeed)
159 write(*,*) 'oobs_pseudoSST: seed for random shuffle of sea-ice points: ', randomSeed
160
161 call oobs_computeObsData(obsData, iceDomainIndexes, iceLons, iceLats, &
162 iceThinning, outputSST, outputFreshWaterST, &
163 outputFileName, datePrint, timePrint, &
164 seaWaterFraction, seaWaterThreshold, inlandWaterPoints)
165 else
166
167 call obs_initialize(obsData, numHeader_max_opt = 0, numBody_max_opt = 0, mpi_local_opt = .true.)
168 call sqlr_writeEmptyPseudoSSTobsFile(obsData, 'SF', outputFileName)
169
170 end if
171
172 if (numberIceCoveredPoints > 0) then
173 deallocate(iceLons)
174 deallocate(iceLats)
175 deallocate(iceDomainIndexes)
176 deallocate(seaWaterFraction)
177 end if
178
179 write(*,*) 'oobs_pseudoSST: done'
180
181 end subroutine oobs_pseudoSST
182
183 !--------------------------------------------------------------------------
184 ! oobs_computeObsData
185 !--------------------------------------------------------------------------
186
187 subroutine oobs_computeObsData(obsData, iceDomainIndexes, iceLons, iceLats, iceThinning, &
188 outputSST, outputFreshWaterST, outputFileName, &
189 datePrint, timePrint, &
190 seaWaterFraction, seaWaterThreshold, inlandWaterPoints)
191 !
192 !:Purpose: pseudo SST data are put into obsSpaceData
193 ! and written into an SQLite file
194
195 implicit none
196
197 ! Arguments:
198 type(struct_obs) , intent(inout) :: obsData ! obsSpaceData
199 integer , intent(in) :: iceDomainIndexes(:)! array of the ice-covered point indexes
200 real(4) , intent(in) :: iceLons(:) ! longitudes of sea ice
201 real(4) , intent(in) :: iceLats(:) ! latitudes of sea ice
202 integer , intent(in) :: iceThinning ! generate pseudo obs in every 'iceThinning' points
203 real(4) , intent(in) :: outputSST ! output SST value for pseudo observations
204 real(4) , intent(in) :: outputFreshWaterST ! output fresh water surface temperature for pseudo obs
205 character(len=*) , intent(in) :: outputFileName
206 integer , intent(in) :: datePrint
207 integer , intent(in) :: timePrint
208 real(4) , intent(in) :: seaWaterFraction(:)! sea water fraction data: 0: fresh water; 1: sea water
209 real(4) , intent(in) :: seaWaterThreshold ! to distinguish inland water from sea water
210 integer , intent(in) :: inlandWaterPoints ! number of inland water points
211
212 ! Locals:
213 real(pre_obsReal) :: obsLon, obsLat
214 real(4) :: obsValue
215 integer :: iceIndex, iceDomainDimension, pseudoObsDimension
216 integer :: codeType, headerIndex
217 integer :: coordinatesIndex, counterThinning, checkInlandWatersCount, checkSeaWatersCount
218 character(len=*), parameter :: myName = 'oobs_computeObsData'
219
220 iceDomainDimension = size(iceDomainIndexes)
221 pseudoObsDimension = floor(real((iceDomainDimension - inlandWaterPoints) / iceThinning)) + inlandWaterPoints
222
223 write(*,*) myName//': sea-ice domain dimension: ', iceDomainDimension
224 write(*,*) myName//': pseudo obs vector dimension: ', pseudoObsDimension
225 write(*,*) myName//': pseudo SST obs will be generated in every ', iceThinning, &
226 ' points of the sea-ice field for sea water, '
227 write(*,*) myName//': and in every point for inland waters, where sea water fraction <= ', seaWaterThreshold
228 write(*,*) myName//': number of inland waters points: ', inlandWaterPoints
229
230 call obs_initialize(obsData, numHeader_max_opt = pseudoObsDimension, &
231 numBody_max_opt = pseudoObsDimension, mpi_local_opt = .true.)
232 codeType = codtyp_get_codtyp('pseudosfc')
233
234 headerIndex = 1
235 counterThinning = iceThinning
236 checkInlandWatersCount = 0
237 checkSeaWatersCount = 0
238
239 do iceIndex = 1, iceDomainDimension
240
241 if (headerIndex > pseudoObsDimension ) cycle
242
243 coordinatesIndex = iceDomainIndexes(iceIndex)
244 obsLon = iceLons(coordinatesIndex)
245 obsLat = iceLats(coordinatesIndex)
246
247 if (seaWaterFraction(coordinatesIndex) <= seaWaterThreshold) then
248 obsValue = (1.0d0 - seaWaterFraction(coordinatesIndex)) * outputFreshWaterST + &
249 seaWaterFraction(coordinatesIndex)* outputSST
250 checkInlandWatersCount = checkInlandWatersCount + 1
251 else
252 if (counterThinning == iceThinning) then
253 obsValue = outputSST
254 checkSeaWatersCount = checkSeaWatersCount + 1
255 counterThinning = 1
256 else
257 counterThinning = counterThinning + 1
258 cycle
259 end if
260 end if
261
262 call obs_setFamily(obsData, 'SF' , headerIndex)
263 call obs_headSet_i(obsData, OBS_ONM, headerIndex, headerIndex)
264 call obs_headSet_i(obsData, OBS_ITY, headerIndex, codeType)
265 call obs_headSet_r(obsData, OBS_LAT, headerIndex, obsLat)
266 call obs_headSet_r(obsData, OBS_LON, headerIndex, obsLon)
267 call obs_bodySet_r(obsData, OBS_VAR, headerIndex, obsValue)
268 call obs_bodySet_i(obsData, OBS_VNM, headerIndex, bufr_sst)
269 call obs_set_c(obsData, 'STID' , headerIndex, 'ABOG')
270 call obs_headSet_i(obsData, OBS_NLV, headerIndex, 1)
271 call obs_headSet_i(obsData, OBS_RLN, headerIndex, headerIndex)
272 call obs_headSet_i(obsData, OBS_DAT, headerIndex, datePrint / 100)
273 call obs_headSet_i(obsData, OBS_ETM, headerIndex, timePrint)
274 headerIndex = headerIndex + 1
275 end do
276
277 call sqlr_writePseudoSSTobs(obsData, 'SF', outputFileName)
278
279 ! Deallocate obsSpaceData
280 call obs_finalize(obsData)
281
282 end subroutine oobs_computeObsData
283
284end module oceanObservations_mod
