1program midas_sstBias
2 !
3 !:Purpose: Main program to compute Sea Surface Temperature (SST)
4 ! satellite observations bias estimate.
5 !
6 ! ---
7 !
8 !:Algorithm: The bias estimation of SST satellite observations is computed
9 ! with respect to insitu observations that are considered unbiased.
10 ! The bias estimation is produced for each sensor separately
11 ! for day and night time.
12 !
13 ! --
14 !
15 ! First, each dataset is put on a regular grid using a small
16 ! search radius (of ~25 km). It is currently a 1800x900 Gaussian grid.
17 ! Second, the bias estimation at every gridpoint is computed as
18 ! an average difference between satellite and insitu observations
19 ! between all collocated valid satellite and insitu observations
20 ! within a larger search radius (of ~1500km).
21 !
22 ! --
23 !
24 ! The resulting bias estimation :math:`B_{a}(k)` at point :math:`k`
25 ! is computed as follows:
26 ! :math:`B_{a}(k) = (1 - w(k)) * B_{b}(k) * \beta + w(k) * B_{a}(k)`,
27 ! where :math:`B_{b}(k)` is a background state of the bias estimation
28 ! computed on the previous day,
29 ! :math:`\beta` is a background term for zero bias in unobserved areas,
30 ! and :math:`w(k)` is a weight which is defined as:
31 ! :math:`w(k) = N_{a}(k) / (N_{a}(k) + N_{b})`,
32 ! where :math:`N_{a}(k)` is the number of observations involved in the
33 ! computation of the current bias estimate :math:`B_{a}(k)` at point :math:`k`
34 ! and :math:`N_{b}` is a parameter for corresponding number
35 ! of observations used to compute the background state :math:`B_{b}(k)`.
36 !
37 !=========================================================== ========================================================
38 ! Input and Output Files Description of file
39 !=========================================================== ========================================================
40 ! ``analysisgrid`` In - File containing the grid where the bias is computed
41 ! ``seaice_analysis`` In - File containing ``LG`` and ``VF`` fields
42 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
43 ! ``searchRadius`` In - 'Large' search radius field to compute biases
44 ! ``trlm_01`` In - Background state of the bias estimation
45 ! ``satellite_bias.fst`` Out - Bias estimations
46 ! ``auxOutput.fst`` Out - Auxiliary output (optional):
47 ! number of observations and weight fields
48 !=========================================================== ========================================================
49 !
50 ! --
51 !
52 !:Synopsis: Below is a summary of the ``SSTbias`` program calling sequence:
53 !
54 ! - **Initial setups:**
55 !
56 ! - Setup horizontal and vertical grid objects for "analysis
57 ! grid" from ``analysisgrid``.
58 !
59 ! - Setup ``obsSpaceData`` object and read observations from
60 ! files: ``inn_setupObs``.
61 !
62 ! - Setup ``columnData`` and ``gridStateVector`` modules.
63 !
64 ! - **Computation**
65 !
66 ! - ``ocm_readMaskFromFile`` get the land-ocean mask
67 !
68 ! - ``oobs_computeObsData`` compute pseudo observation values
69 ! and their coordinates and save them in SQLite files.
70 !
71 ! - ``sstb_getGriddedObs`` get all datasets on a regular grid
72 !
73 ! - ``sstb_getGriddedBias`` compute bias estimation for each sensor,
74 ! for day time or night time on a regular grid,
75 ! and save the results into an output standard file.
76 ! --
77 !
78 !:Options: `List of namelist blocks <../namelists_in_each_program.html#sstbias>`_
79 ! that can affect the ``SSTbias`` program.
80 !
81 ! * The use of ``SSTbias`` program is controlled by the namelist block
82 ! ``&namSSTbiasEstimate`` read by the ``SSTbias`` program.
83 !
84 ! * ``iceFractionThreshold`` the sea-ice fraction threshold to define
85 ! the presence of ice
86 !
87 ! * ``searchRadius`` horizontal search radius for observation gridding
88 !
89 ! * ``maxBias`` max allowed insitu-satellite difference in degrees
90 !
91 ! * ``numberPointsBG`` :math:`N_{b}`, number of points to compute
92 ! the background bias estimation
93 !
94 ! * ``sensorList`` name of sensor
95 !
96 ! * ``weightMin`` minimum value of weight
97 !
98 ! * ``weightMax`` maximum value of weight
99 !
100 ! * ``saveAuxFields`` to store or not auxiliary fields: nobs and weight
101 !
102 ! * ``bgTermZeroBias`` background term to zero bias
103 !
104 ! --
105 !
106 use version_mod
107 use ramDisk_mod
108 use utilities_mod
109 use midasMpi_mod
110 use mathPhysConstants_mod
111 use horizontalCoord_mod
112 use verticalCoord_mod
113 use timeCoord_mod
114 use obsSpaceData_mod
115 use gridStateVector_mod
116 use obsFiles_mod
117 use innovation_mod
118 use sstBias_mod
119 use columnData_mod
120
121 implicit none
122
123 integer, external :: exdb, exfin, fnom, fclos, get_max_rss
124 integer :: ierr, istamp
125
126 type(struct_obs), target :: obsSpaceData
127 type(struct_hco), pointer :: hco_anl => null()
128 type(struct_vco), pointer :: vco_anl => null()
129 character(len=48),parameter :: obsMpiStrategy = 'LIKESPLITFILES'
130 character(len=48),parameter :: varMode = 'analysis'
131 type(struct_columnData) :: column ! column data
132 integer :: dateStampFromObs
133
134 istamp = exdb('SSTBIASESTIMATION','DEBUT','NON')
135
136 call ver_printNameAndVersion('SSTbias','SST Bias Estimation')
137
138 ! MPI initialization
139 call mmpi_initialize
140
141 call tmg_init(mmpi_myid, 'TMG_INFO')
142
143 call utl_tmg_start(0,'Main')
144
145 ! 1. Top level setup
146
147 call ram_setup()
148
149 ! Do initial set up
150 call SSTbias_setup('VAR') ! obsColumnMode
151
152 call sstb_computeBias(obsSpaceData, hco_anl, vco_anl)
153
154 ! Deallocate copied obsSpaceData
155 call obs_finalize(obsSpaceData)
156 call col_deallocate(column)
157
158 ! 3. Job termination
159
160 istamp = exfin('SSTBIAS','FIN','NON')
161
162 call utl_tmg_stop(0)
163
164 call tmg_terminate(mmpi_myid, 'TMG_INFO')
165
166 call rpn_comm_finalize(ierr)
167
168 contains
169
170 subroutine SSTbias_setup(obsColumnMode)
171 !
172 ! :Purpose: Control of the preprocessing of bias estimation
173 !
174 implicit none
175
176 ! Arguments:
177 character(len=*), intent(in) :: obsColumnMode
178
179 ! Locals:
180 character(len=*), parameter :: gridFile = './analysisgrid'
181
182 write(*,*) ''
183 write(*,*) '-------------------------------------------------'
184 write(*,*) '-- Starting subroutine SSTbias_setup --'
185 write(*,*) '-------------------------------------------------'
186
187 !
188 !- Initialize the Temporal grid and set dateStamp from env variable
189 !
190 call tim_setup()
191
192 !
193 !- Initialize observation file names and set dateStamp
194 !
195 call obsf_setup(dateStampFromObs, varMode)
196 if (tim_getDateStamp() == 0) then
197 if (dateStampFromObs > 0) then
198 ! use dateStamp from obs if not set by env variable
199 call tim_setDateStamp(dateStampFromObs)
200 else
201 call utl_abort('SSTbias_setup: DateStamp was not set')
202 end if
203 end if
204
205 !
206 !- Initialize constants
207 !
208 if(mmpi_myid == 0) call mpc_printConstants(6)
209
210 !
211 !- Initialize variables of the model states
212 !
213 call gsv_setup
214 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
215
216 !
217 !- Initialize the Analysis grid
218 !
219 if(mmpi_myid == 0) write(*,*)''
220 if(mmpi_myid == 0) write(*,*) 'SSTbias_setup: Set hco parameters for analysis grid'
221 call hco_SetupFromFile(hco_anl, gridFile, 'GRID' ) ! IN
222
223 !
224 !- Initialisation of the analysis grid vertical coordinate from analysisgrid file
225 !
226 call vco_SetupFromFile(vco_anl, & ! OUT
227 gridFile ) ! IN
228
229 call col_setVco(column, vco_anl)
230 !
231 !- Setup and read observations
232 !
233 call inn_setupObs(obsSpaceData, hco_anl, obsColumnMode, obsMpiStrategy, varMode) ! IN
234 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
235
236 !- Basic setup of columnData module
237 call col_setup
238
239 !- Memory allocation for background column data
240 call col_allocate(column, obs_numHeader(obsSpaceData), mpiLocal_opt = .true.)
241
242 if(mmpi_myid == 0) write(*,*) 'SSTbias_setup: done.'
243
244 end subroutine SSTbias_setup
245
246end program midas_sstBias
