1program midas_pseudoSSTobs
2 !
3 !:Purpose: Main program to produce pseudo SST observations
4 ! in ice-covered areas. Pseudo SST observations are needed
5 ! to prevent the propagation of analysis increments to
6 ! the ice-covered areas, that may result in undesirable sea-ice melting.
7 !
8 ! ---
9 !
10 !:Algorithm: Pseudo SST observations are assigned to the ice-covered
11 ! water points.
12 ! First, a global sea-ice analysis is read.
13 ! The sea-ice analysis file contains a mandatory sea-water
14 ! fraction field.
15 ! The grid and land-ocean mask are read
16 ! from the ``analysisgrid`` file.
17 !
18 ! --
19 !
20 ! Second, the number of ice-covered water points, including
21 ! concerned inland water points, are computed.
22 ! If the number of ice-covered water points is zero,
23 ! an empty observation SQLite file is created.
24 ! If not, the computation of pseudo observations starts.
25 !
26 ! --
27 !
28 ! First, the index array of ice-covered water points are
29 ! randomly shuffled to prevent the insertion of pseudo
30 ! observations at the same locations
31 ! that would lead to spatial correlation of observations.
32 ! Second, the pseudo observations of sea surface temperature :math:`T`
33 ! are inserted at every ice-covered inland water point :math:`k`,
34 ! where the value of observations is computed as follows:
35 ! :math:`T(k)=(1 - w(k)) * T_{fw} + w(k) * T_{s}`,
36 ! where :math:`w(k)` is the sea-water fraction at the point :math:`k`,
37 ! :math:`T_{fw}` is the temperature of fresh water below the ice,
38 ! :math:`T_{s}` is a temperature of the sea water below the ice.
39 ! The pseudo observations are inserted into every :math:`N`-th point
40 ! of sea water ice-covered points,
41 ! where the value of observation is defined as
42 ! :math:`T_{s}`.
43 !
44 ! --
45 !
46 ! The computed observation values along with the corresponding
47 ! coordinates are put into ``obsSpaceData``.
48 ! Finally, output SQLite files are created.
49 !
50 ! --
51 !
52 !=========================================================== ======================================================
53 ! Input and Output Files Description of file
54 !=========================================================== ======================================================
55 ! ``analysisgrid`` In - File defining sea-ice global grid
56 ! ``seaice_analysis`` In - File containing ``LG`` and ``VF`` fields
57 ! ``obsfiles_sst_pseudo.updated/obssst_pseudo_$NNNN_$NNNN`` Out - Pseudo obs file for each MPI task
58 !=========================================================== ======================================================
59 !
60 ! --
61 !
62 !:Synopsis: Below is a summary of the ``pseudoSSTobs`` program calling sequence:
63 !
64 ! - **Initial setups:**
65 !
66 ! - Setup horizontal and vertical grid objects for "analysis
67 ! grid" from ``analysisgrid``.
68 !
69 ! - Setup ``obsSpaceData`` object.
70 !
71 ! - Setup ``gridStateVector`` module.
72 !
73 ! - **Computation**
74 !
75 ! - ``utl_randomOrderInt`` random shuffle the ice-covered point indices
76 !
77 ! - ``oobs_computeObsData`` compute pseudo observation locations and values
78 ! and save them in SQLite files.
79 !
80 ! --
81 !
82 !:Options: `List of namelist blocks <../namelists_in_each_program.html#pseudoSSTobs>`_
83 ! that can affect the ``pseudoSSTobs`` program.
84 !
85 ! * The use of ``pseudoSSTobs`` program is controlled by the namelist block
86 ! ``&pseudoSSTobs`` read by the ``pseudoSSTobs`` program.
87 !
88 ! * ``iceFractionThreshold`` the sea-ice fraction threshold to define
89 ! the presence of ice at each particular point
90 !
91 ! * ``outputSST`` the value of :math:`T_{s}` in K;
92 !
93 ! * ``outputFreshWaterST`` the value of :math:`T_{fw}` in K;
94 !
95 ! * ``seaiceThinning`` pseudo observations are inserted into each :math:`N`-th point,
96 ! this parameter controls the observation thinning
97 !
98 ! * ``outputFileName`` controls the output file names
99 !
100 ! * ``etiket`` etiket to put into the table "resume" of output SQLite file
101 !
102 ! * ``seaWaterThreshold`` a threshold to distinguish sea and fresh water
103 !
104 ! --
105 !
106 !========================== ================ ==============================================================
107 ! Module Namelist Description of what is controlled
108 !========================== ================ ==============================================================
109 ! ``oceanObservations_mod`` ``pseudoSSTobs`` parameters of pseudo SST observations
110 !========================== ================ ==============================================================
111 !
112 use version_mod
113 use ramDisk_mod
114 use utilities_mod
115 use horizontalCoord_mod
116 use verticalCoord_mod
117 use midasMpi_mod
118 use oceanObservations_mod
119 use gridStateVector_mod
120 use obsSpaceData_mod
121
122 implicit none
123
124 integer, external :: exdb, exfin, fnom, fclos, get_max_rss
125 integer :: ierr, istamp, nulnam
126
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
132 ! namelist variables
133 real(4) :: iceFractionThreshold ! consider no ice condition below this threshold
134 real(4) :: outputSST ! output SST value for pseudo observations
135 real(4) :: outputFreshWaterST ! output fresh water surface temperature for pseudo obs.
136 integer :: seaiceThinning ! generate pseudo obs in every 'seaiceThinning' points
137 character(len=100) :: outputFileName ! name of the file containing the generated observations
138 real(4) :: seaWaterThreshold ! to distinguish inland water from sea water
139
140 namelist /pseudoSSTobs/ iceFractionThreshold, outputSST, outputFreshWaterST, seaiceThinning, &
141 outputFileName, seaWaterThreshold
142
143 istamp = exdb('pseudoSSTobs','DEBUT','NON')
144
145 call ver_printNameAndVersion('pseudoSSTobs','Generation of pseudo SST observations')
146
147 ! MPI initialization
148 call mmpi_initialize
149
150 call tmg_init(mmpi_myid, 'TMG_INFO')
151
152 call utl_tmg_start(0,'Main')
153
154 ! 1. Top level setup
155
156 call ram_setup()
157
158 ! Do initial set up
159 call pseudoSSTobs_setup()
160
161 call oobs_pseudoSST(hco_anl, vco_anl, iceFractionThreshold, outputSST, outputFreshWaterST, &
162 seaiceThinning, outputFileName, seaWaterThreshold)
163
164 ! 3. Job termination
165
166 istamp = exfin('pseudoSSTobs','FIN','NON')
167
168 call utl_tmg_stop(0)
169
170 call tmg_terminate(mmpi_myid, 'TMG_INFO')
171
172 call rpn_comm_finalize(ierr)
173
174 contains
175
176 subroutine pseudoSSTobs_setup()
177 !
178 ! :Purpose: Control of the preprocessing of pseudo SST obs
179 !
180 implicit none
181
182 ! Locals:
183 character(len=*), parameter :: myName = 'pseudoSSTobs_setup'
184 character(len=*), parameter :: gridFile = './analysisgrid'
185
186 write(*,*) ''
187 write(*,*) '-------------------------------------------------'
188 write(*,*) '-- Starting subroutine '//myName//' --'
189 write(*,*) '-------------------------------------------------'
190
191 ! Setting default namelist variable values
192 iceFractionThreshold = 0.05
193 outputSST = 271.4
194 outputFreshWaterST = 271.4
195 outputFileName = ''
196 seaiceThinning = 5
197 seaWaterThreshold = 0.0
198
199 ! Read the namelist
200 nulnam = 0
201 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
202 read(nulnam, nml = pseudoSSTobs, iostat = ierr)
203 if (ierr /= 0) call utl_abort(myName//': Error reading namelist')
204 if (mmpi_myid == 0) write(*, nml = pseudoSSTobs)
205 ierr = fclos(nulnam)
206
207 write(*,*)''
208 write(*,*) myName//': output SST globally: ', outputSST
209 write(*,*) myName//': output fresh water surface temperature globally: ', outputFreshWaterST
210 write(*,*) myName//': sea-ice fraction threshold: ', iceFractionThreshold
211 write(*,*) myName//': sea water fraction threshold: ', seaWaterThreshold
212 write(*,*) myName//': pseudo SST obs will be generated in every ', seaiceThinning, ' points of the sea-ice field'
213 write(*,*) myName//': output file name: ', outputFileName
214 !
215 !- Initialize the Analysis grid
216 !
217 if(mmpi_myid == 0) write(*,*)''
218 if(mmpi_myid == 0) write(*,*) myName//': Set hco parameters for analysis grid'
219 call hco_SetupFromFile(hco_anl, gridFile, 'ANALYSIS') ! IN
220
221 !
222 !- Initialisation of the analysis grid vertical coordinate from analysisgrid file
223 !
224 call vco_SetupFromFile(vco_anl, & ! OUT
225 gridFile) ! IN
226
227 !- Initialize variables of the model states
228 !
229 call gsv_setup
230
231 call obs_class_initialize('VAR')
232
233 if(mmpi_myid == 0) write(*,*) myName//': done.'
234
235 end subroutine pseudoSSTobs_setup
236
237end program midas_pseudoSSTobs
