1module fsoi_mod
2 ! MODULE fsoi_mod (prefix='fso' category='1. High-level functionality')
3 !
4 !:Purpose: Observation impact (FSOI) library.
5 !
6 use midasMpi_mod
7 use codePrecision_mod
8 use horizontalCoord_mod
9 use verticalCoord_mod
10 use obsSpaceData_mod
11 use gridStateVector_mod
12 use gridStateVectorFileIO_mod
13 use columnData_mod
14 use controlVector_mod
15 use bmatrix_mod
16 use bmatrixensemble_mod
17 use stateToColumn_mod
18 use obsOperators_mod
19 use rMatrix_mod
20 use MathPhysConstants_mod
21 use quasinewton_mod
22 use costFunction_mod
23 use tovsNL_mod
24 use timeCoord_mod
25 use utilities_mod
26 use calcHeightAndPressure_mod
27 use rttov_const, only: inst_name, platform_name
28 implicit none
29 save
30 private
31
32 ! public subroutines and functions
33 public :: fso_setup, fso_ensemble
34
35 ! module private variables
36 type(struct_obs), pointer :: obsSpaceData_ptr
37 type(struct_columnData), pointer :: columnTrlOnAnlIncLev_ptr
38 type(struct_columnData), pointer :: column_ptr
39 real(8),allocatable :: vhat(:)
40 integer,external :: get_max_rss
41 integer :: fso_nsim, nvadim_mpilocal
42
43 type(struct_hco), pointer :: hco_anl => null()
44
45 ! namelist variables
46 integer :: nvamaj ! number of vector pairs to store in memory for Hessian approximation
47 integer :: nitermax ! maximum number of minimization iterations
48 integer :: nsimmax ! maximum number of cost function evaluations during minimization
49 real(8) :: leadTime ! lead time of forecast (in hours)
50 real(8) :: repsg ! relative gradient amplitude used as stopping criteria
51 real(8) :: rdf1fac ! factor applied to initial cost function value to approximate final value
52 real(8) :: latMinNorm ! minimum latitude for area included in forecast error norm (in degrees)
53 real(8) :: latMaxNorm ! maximum latitude for area included in forecast error norm (in degrees)
54 real(8) :: lonMinNorm ! minimum longitude for area included in forecast error norm (in degrees)
55 real(8) :: lonMaxNorm ! maximum longitude for area included in forecast error norm (in degrees)
56 logical :: includeUVnorm ! choose to include winds in forecast error norm
57 logical :: includeTTnorm ! choose to include temperature in forecast error norm
58 logical :: includeP0norm ! choose to include surface pressure in forecast error norm
59 logical :: includeHUnorm ! choose to include humidity in forecast error norm
60 logical :: includeTGnorm ! choose to include surface skin temperature in forecast error norm
61 character(len=256) :: forecastPath ! relative path where forecast files are stored
62 character(len=4) :: fsoMode ! type of FSOI algorithm: can be 'HFSO' or 'EFSO'
63 logical :: StratoNorm ! choose for forecast error norm from 100hPa to 1hPa,default from surface to 100hPa
64
65 contains
66
67 !--------------------------------------------------------------------------
68 ! fso_setup
69 !--------------------------------------------------------------------------
70 subroutine fso_setup(hco_anl_in)
71 !
72 !:Purpose: Initialise the FSOI module: read the namelist and initialise
73 ! global variables and structure
74 !
75 implicit none
76
77 ! Arguments:
78 type(struct_hco), pointer, intent(in) :: hco_anl_in
79
80 ! Locals:
81 integer :: ierr,nulnam
82 integer :: fnom,fclos
83
84 NAMELIST /NAMFSO/leadTime, nvamaj, nitermax, nsimmax
85 NAMELIST /NAMFSO/repsg, rdf1fac, forecastPath, fsoMode
86 NAMELIST /NAMFSO/latMinNorm, latMaxNorm, lonMinNorm, lonMaxNorm
87 NAMELIST /NAMFSO/includeUVnorm, includeTTnorm, includeP0norm, includeHUnorm, includeTGnorm
88 NAMELIST /NAMFSO/StratoNorm
89
90 ! set default values for namelist variables
91 leadtime = 12.0d0
92 nvamaj = 6
93 nitermax = 100
94 nsimmax = 120
95 repsg = 1d-5
96 rdf1fac = 0.25d0
97 latMinNorm = -95.0d0
98 latMaxNorm = 95.0d0
99 lonMinNorm = -185.0d0
100 lonMaxNorm = 365.0d0
101 includeUVnorm=.true.
102 includeTTnorm=.true.
103 includeP0norm=.true.
104 includeHUnorm=.false.
105 includeTGnorm=.false.
106 forecastPath = './forecasts'
107 fsoMode = 'HFSO'
108 StratoNorm = .false.
109
110 ! read in the namelist NAMFSO
111 nulnam = 0
112 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
113 read(nulnam,nml=namfso,iostat=ierr)
114 if(ierr /= 0) call utl_abort('fso_setup: Error reading namelist')
115 write(*,nml=namfso)
116 ierr = fclos(nulnam)
117
118 if (trim(fsoMode) /= 'HFSO' .and. trim(fsoMode) /= 'EFSO') then
119 call utl_abort('fso_setup: Invalid value of fsoMode. Must be HFSO or EFSO')
120 end if
121
122 ! convert Latmin,max and Lonmin,max from degres into RAD
123 latMinNorm = latMinNorm*MPC_RADIANS_PER_DEGREE_R8
124 latMaxNorm = latMaxNorm*MPC_RADIANS_PER_DEGREE_R8
125 lonMinNorm = lonMinNorm*MPC_RADIANS_PER_DEGREE_R8
126 lonMaxNorm = lonMaxNorm*MPC_RADIANS_PER_DEGREE_R8
127
128 call ben_setFsoLeadTime(leadTime)
129 fso_nsim = 0
130
131 hco_anl => hco_anl_in
132
133 end subroutine fso_setup
134
135 !--------------------------------------------------------------------------
136 ! fso_ensemble
137 !--------------------------------------------------------------------------
138 subroutine fso_ensemble(columnTrlOnAnlIncLev,obsSpaceData)
139 !
140 !:Purpose: Perform forecast sensitivity to observation calculation using
141 ! ensemble approach
142 !
143 implicit none
144
145 ! Arguments:
146 type(struct_columnData), target, intent(inout) :: columnTrlOnAnlIncLev
147 type(struct_obs), target, intent(inout) :: obsSpaceData
148
149 ! Locals:
150 type(struct_columnData),target :: column
151 type(struct_gsv) :: statevector_FcstErr, statevector_fso, statevector_HUreference
152 type(struct_vco), pointer :: vco_anl
153 real(8),allocatable :: ahat(:), zhat(:)
154 integer :: dateStamp_fcst, dateStamp
155 integer :: headerIndex, bodyIndexBeg, bodyIndexEnd, bodyIndex
156 real(8) :: fso_ori, fso_fin
157
158 if (mmpi_myid == 0) write(*,*) 'fso_ensemble: starting'
159
160 vco_anl => col_getVco(columnTrlOnAnlIncLev)
161
162 nvadim_mpilocal = cvm_nvadim
163
164 ! initialize column object for storing "increment"
165 call col_setVco(column,col_getVco(columnTrlOnAnlIncLev))
166 call col_allocate(column,col_getNumCol(columnTrlOnAnlIncLev))
167
168 ! compute dateStamp_fcst
169 call incdatr(dateStamp_fcst, tim_getDatestamp(), leadTime)
170 dateStamp = tim_getDatestamp()
171 write(*,*) 'fso_ensemble: analysis datestamp = ',dateStamp
172 write(*,*) 'fso_ensemble: forecast datestamp = ',dateStamp_fcst
173
174 ! allocate control vector related arrays (these are all mpilocal)
175 allocate(ahat(nvadim_mpilocal))
176 allocate(vhat(nvadim_mpilocal))
177 allocate(zhat(nvadim_mpilocal))
178
179 ! initialize control vector related arrays to zero
180 ahat(:) = 0.0d0
181 vhat(:) = 0.0d0
182
183 ! for statevector_FcstErr
184 call gsv_allocate(statevector_FcstErr, 1, hco_anl, vco_anl, &
185 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
186 allocHeight_opt=.false., allocPressure_opt=.false.)
187
188 ! for statevector_fso
189 call gsv_allocate(statevector_fso, tim_nstepobsinc, hco_anl, vco_anl, &
190 dataKind_opt=pre_incrReal, &
191 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true.)
192
193 ! for statevector_HUreference (verifying analysis)
194 call gsv_allocate(statevector_HUreference, 1, hco_anl, vco_anl, &
195 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
196 allocHeight_opt=.false., allocPressure_opt=.false.)
197
198 ! compute forecast error = C * (error_t^fa + error_t^fb)
199 call calcFcstError(columnTrlOnAnlIncLev,statevector_FcstErr,statevector_HUreference)
200
201
202 ! compute vhat = B_t^T/2 * C * (error_t^fa + error_t^fb)
203 call bmat_sqrtBT(vhat, nvadim_mpilocal, statevector_FcstErr, useFSOFcst_opt = .true., &
204 stateVectorRef_opt=statevector_HUreference)
205
206 if (mmpi_myid == 0) write(*,*) 'fso: B_t^T/2 * C * (error_t^fa + error_t^fb) max,min:', &
207 maxval(vhat),minval(vhat)
208
209 if( trim(fsoMode) == 'HFSO' ) then
210 call minimize(nvadim_mpilocal, zhat, column, columnTrlOnAnlIncLev, obsSpaceData)
211 ahat = zhat + vhat
212 call bmat_sqrtB(ahat, nvadim_mpilocal, statevector_fso)
213 elseif( trim(fsoMode) == 'EFSO' ) then
214 call bmat_sqrtB(vhat, nvadim_mpilocal, statevector_fso)
215 end if
216
217 ! Compute yhat = [R^-1 H B^1/2 ahat], and put in OBS_FSO
218 call s2c_tl(statevector_fso,column,columnTrlOnAnlIncLev,obsSpaceData) ! put in column H_horiz B^1/2 ahat
219 call oop_Htl(column,columnTrlOnAnlIncLev,obsSpaceData,1) ! Save as OBS_WORK: H_vert H_horiz B^1/2 vhat = H B^1/2 ahat
220 call rmat_RsqrtInverseAllObs(obsSpaceData,OBS_FSO,OBS_WORK) ! Save as OBS_FSO : R**-1/2 H B^1/2 ahat
221 call rmat_RsqrtInverseAllObs(obsSpaceData,OBS_FSO,OBS_FSO) ! Save as OBS_FSO : R**-1 H B^1/2 ahat\
222
223 do headerIndex = 1, obs_numHeader(obsSpaceData)
224
225 bodyIndexBeg = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
226 bodyIndexEnd = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + bodyIndexBeg - 1
227
228 do bodyIndex = bodyIndexBeg, bodyIndexEnd
229 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
230 fso_ori = obs_bodyElem_r(obsSpaceData,OBS_FSO,bodyIndex)
231 fso_fin = fso_ori * obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)
232 call obs_bodySet_r(obsSpaceData,OBS_FSO,bodyIndex, fso_fin)
233 end if
234 end do
235
236 end do
237
238 ! print out the information of total FSO for each family
239 call sumFSO(obsSpaceData)
240
241 ! deallocate the control vector related arrays
242 deallocate(ahat)
243 deallocate(vhat)
244 deallocate(zhat)
245 call col_deallocate(column)
246
247 if (mmpi_myid == 0) write(*,*) 'fso_ensemble: Finished'
248
249 end subroutine fso_ensemble
250
251 !--------------------------------------------------------------------------
252 ! calcFcstError
253 !--------------------------------------------------------------------------
254 subroutine calcFcstError(columnTrlOnAnlIncLev,statevector_out,statevector_verifAnalysis)
255 !
256 !:Purpose: Reads the forecast from background and analysis, the verifying
257 ! analysis based on these inputs, calculates the Forecast error
258 !
259 implicit none
260
261 ! Arguments:
262 type(struct_columnData), target, intent(in) :: columnTrlOnAnlIncLev
263 type(struct_gsv) , target, intent(inout) :: statevector_out
264 type(struct_gsv) , target, intent(inout) :: statevector_verifAnalysis
265
266 ! Locals:
267 type(struct_gsv) :: statevector_fa, statevector_fb, statevector_a
268 character(len=256) :: fileName_fa, fileName_fb, fileName_a
269 logical :: faExists
270 type(struct_gsv) :: statevector_tempfa, statevector_tempfb
271 type(struct_vco), pointer :: vco_anl
272 integer :: dateStamp_fcst, dateStamp
273
274 vco_anl => col_getVco(columnTrlOnAnlIncLev)
275
276 ! compute dateStamp_fcst
277 call incdatr(dateStamp_fcst, tim_getDatestamp(), leadTime)
278 dateStamp = tim_getDatestamp()
279 write(*,*) 'fso_ensemble: analysis datestamp = ',dateStamp
280 write(*,*) 'fso_ensemble: forecast datestamp = ',dateStamp_fcst
281
282 ! read forecasts from the analysis and background state
283 fileName_fa = trim(forecastPath) // '/forecast_a'
284 inquire(file=trim(fileName_fa),exist=faExists)
285 write(*,*) 'faExists', faExists
286 call gsv_allocate(statevector_fa, 1, hco_anl, vco_anl, &
287 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
288 hInterpolateDegree_opt='LINEAR', &
289 allocHeight_opt=.false., allocPressure_opt=.false.)
290 call gio_readFromFile(statevector_fa, fileName_fa, ' ', 'P', containsFullField_opt=.true.)
291
292 !for statevecotr_tempfa
293 call gsv_allocate(statevector_tempfa, 1, hco_anl, vco_anl, &
294 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
295 allocHeight_opt=.false., allocPressure_opt=.false.)
296
297 !for statevector_fb
298 fileName_fb = trim(forecastPath) // '/forecast_b'
299 call gsv_allocate(statevector_fb, 1, hco_anl, vco_anl, &
300 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
301 hInterpolateDegree_opt='LINEAR', &
302 allocHeight_opt=.false., allocPressure_opt=.false.)
303 call gio_readFromFile(statevector_fb, fileName_fb, ' ', 'P', containsFullField_opt=.true.)
304
305 !for statevecotr_tempfb
306 call gsv_allocate(statevector_tempfb, 1, hco_anl, vco_anl, &
307 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
308 allocHeight_opt=.false., allocPressure_opt=.false.)
309
310 ! read verifying analysis
311 fileName_a = trim(forecastPath) // '/analysis'
312 call gsv_allocate(statevector_a, 1,hco_anl, vco_anl, &
313 datestamp_opt=datestamp_fcst, mpi_local_opt=.true., &
314 hInterpolateDegree_opt='LINEAR', &
315 allocHeight_opt=.false., allocPressure_opt=.false.)
316 call gio_readFromFile(statevector_a, fileName_a, ' ', 'A', containsFullField_opt=.true.)
317
318 ! compute error of both forecasts (overwrite forecasts with error)
319 call gsv_add(statevector_a, statevector_fa, -1.0d0)
320 call gsv_add(statevector_a, statevector_fb, -1.0d0)
321
322 call gsv_copy(statevector_fa,statevector_tempfa)
323 call gsv_copy(statevector_fb,statevector_tempfb)
324 call multEnergyNorm(statevector_tempfa, statevector_a, &
325 latMinNorm, latMaxNorm, lonMinNorm, lonMaxNorm, &
326 includeUVnorm, includeTTnorm, includeP0norm, &
327 includeHUnorm, includeTGnorm, StratoNorm) ! use analysis as reference state
328 call multEnergyNorm(statevector_tempfb, statevector_a, &
329 latMinNorm, latMaxNorm, lonMinNorm, lonMaxNorm, &
330 includeUVnorm, includeTTnorm, includeP0norm, &
331 includeHUnorm, includeTGnorm, StratoNorm) ! use analysis as reference state
332
333 ! compute error Norm = C * (error_t^fa + error_t^fb)
334 call gsv_add(statevector_fa, statevector_fb, 1.0d0)
335 call multEnergyNorm(statevector_fb, statevector_a, &
336 latMinNorm, latMaxNorm, lonMinNorm, lonMaxNorm, &
337 includeUVnorm, includeTTnorm, includeP0norm, &
338 includeHUnorm, includeTGnorm, StratoNorm) ! use analysis as reference state
339 call gsv_copy(statevector_fb,statevector_out)
340 call gsv_copy(statevector_a,statevector_verifAnalysis)
341
342 end subroutine calcFcstError
343
344 !--------------------------------------------------------------------------
345 ! minimize
346 !--------------------------------------------------------------------------
347 subroutine minimize(nvadim,zhat,column,columnTrlOnAnlIncLev,obsSpaceData)
348 !
349 !:Purpose: Performs HFSO quasi-Newton minimization
350 !
351 implicit none
352
353 ! Arguments:
354 integer, intent(in) :: nvadim
355 real(8), dimension(nvadim), intent(out) :: zhat
356 type(struct_columnData), target, intent(in) :: column
357 type(struct_columnData), target, intent(in) :: columnTrlOnAnlIncLev
358 type(struct_obs), target, intent(in) :: obsSpaceData
359
360 ! Locals:
361 integer :: nulout = 6
362 integer :: imode, itermax, isimmax, indic, nmtra
363 integer :: impres, iztrl(10), intunused(1)
364 real :: rspunused(1)
365 real(8) :: zjsp, zxmin, zdf1, zeps, dlgnorm, dlxnorm,zspunused(1)
366 real(8),allocatable :: gradJ(:), vatra(:)
367
368 call utl_tmg_start(90,'--Minimization')
369
370 if (mmpi_myid == 0) write(*,*) 'minimize: starting'
371
372 nmtra = (4 + 2*nvamaj)*nvadim
373 write(*,'(4X,"NVAMAJ = ",I3,/5X,"NMTRA =",I14)') nvamaj,nmtra
374
375 columnTrlOnAnlIncLev_ptr => columnTrlOnAnlIncLev
376 column_ptr => column
377 obsSpaceData_ptr => obsSpaceData
378
379 allocate(gradJ(nvadim))
380 allocate(vatra(nmtra))
381
382 gradJ(:) = 0.0d0
383 vatra(:) = 0.0d0
384 zhat(:) = 0.0d0
385
386 ! Compute zhat by performing variational minimization
387 ! Set-up for the minimization
388 if (mmpi_myid == 0) then
389 impres = 5
390 else
391 impres = 0
392 end if
393
394 imode = 0
395 zeps = repsg
396 itermax = nitermax
397 isimmax = nsimmax
398 zxmin = epsilon(zxmin)
399 ! initial gradient calculation
400 indic = 2
401 call utl_tmg_start(91,'----QuasiNewton')
402 call simvar(indic,nvadim,zhat,zjsp,gradJ)
403 call utl_tmg_stop(91)
404 zdf1 = rdf1fac * abs(zjsp)
405
406 ! print amplitude of initial gradient and cost function value
407 call prscal(nvadim,gradJ,gradJ,dlgnorm)
408 dlgnorm = dsqrt(dlgnorm)
409 call prscal(nvadim,zhat,zhat,dlxnorm)
410 dlxnorm = dsqrt(dlxnorm)
411 write(*,*)' |X| = ', dlxnorm
412 write(*,'(/4X,"J(X) = ",G23.16,4X,"|Grad J(X)| = ",G23.16)') zjsp, dlgnorm
413 write(*,'(//,10X," Minimization QNA_N1QN3 starts ...",/ &
414 10x,"DXMIN =",G23.16,2X,"DF1 =",G23.16,2X,"EPSG =",G23.16 &
415 /,10X,"IMPRES =",I3,2X,"NITER = ",I3,2X,"NSIM = ",I3)') zxmin,zdf1,zeps,impres,itermax,isimmax
416
417 ! Do the minimization
418 call utl_tmg_start(91,'----QuasiNewton')
419 call qna_n1qn3(simvar, prscal, dcanonb, dcanab, nvadim, zhat, &
420 zjsp, gradJ, zxmin, zdf1, zeps, impres, nulout, imode, &
421 itermax,isimmax, iztrl, vatra, nmtra, intunused, rspunused, &
422 zspunused)
423 call utl_tmg_stop(91)
424 call fool_optimizer(obsSpaceData)
425
426 write(*,'(//,20X,20("*"),2X &
427 ,/,20X," Minimization ended with MODE:",I4 &
428 ,/,20X," Total number of iterations:",I4 &
429 ,/,20X," Total number of simulations:",I4)' ) imode,itermax,isimmax
430
431 if (mmpi_myid == 0) write(*,*) 'minimize: Finished'
432
433 deallocate(vatra)
434 deallocate(gradJ)
435
436 call utl_tmg_stop(90)
437
438 end subroutine minimize
439
440 !--------------------------------------------------------------------------
441 ! sumFSO
442 !--------------------------------------------------------------------------
443 subroutine sumFSO(obsSpaceData)
444 !
445 !:Purpose: Print out the information of total FSO for each family
446 !
447 implicit none
448
449 ! Arguments:
450 type(struct_obs), intent(in) :: obsSpaceData
451
452 ! Locals:
453 real(8) :: pfso_1
454 integer :: bodyIndex,itvs,isens,headerIndex
455 integer :: bodyIndexBeg, bodyIndexEnd
456 integer, parameter :: numFamily = 10
457 character(len=2), parameter :: familyList(numFamily) = (/'UA','AI','SF','SC','TO','SW','PR','RO','GP','CH'/)
458 real(8) :: tfso(numFamily), tfsotov_sensors(tvs_nsensors),totFSO
459 integer :: numAss_local(numFamily), numAss_global(numFamily)
460 integer :: numAss_sensors_loc(tvs_nsensors), numAss_sensors_glb(tvs_nsensors)
461 integer :: ierr, familyIndex
462
463 if (mmpi_myid == 0) write(*,*) 'sumFSO: Starting'
464
465 ! initialize
466 do familyIndex = 1, numFamily
467 tfso(familyIndex) = 0.d0
468 numAss_local(familyIndex) = 0
469 numAss_global(familyIndex) = 0
470 end do
471
472 tfsotov_sensors(:) = 0.d0
473 numAss_sensors_loc(:) = 0
474 numAss_sensors_glb(:) = 0
475 totFSO = 0.d0
476
477 do familyIndex = 1, numFamily
478 do bodyIndex = 1, obs_numbody(obsSpaceData)
479
480 pfso_1 = obs_bodyElem_r(obsSpaceData,OBS_FSO,bodyIndex)
481 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == 1 ) then
482 ! FSO for each family
483 if (obs_getFamily(obsSpaceData,bodyIndex_opt=bodyIndex) == familyList(familyIndex) ) then
484 tfso(familyIndex) = tfso(familyIndex) + pfso_1
485 numAss_local(familyIndex) = numAss_local(familyIndex) + 1
486 end if
487 end if
488
489 end do
490 end do
491
492 do itvs = 1, tvs_nobtov
493 headerIndex = tvs_headerIndex(itvs)
494 if (headerIndex > 0 ) then
495 bodyIndexBeg = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
496 bodyIndexEnd = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + bodyIndexBeg - 1
497 do bodyIndex = bodyIndexBeg, bodyIndexEnd
498 pfso_1 = obs_bodyElem_r(obsSpaceData,OBS_FSO,bodyIndex)
499 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == 1 ) then
500 isens = tvs_lsensor (itvs)
501 tfsotov_sensors(isens) = tfsotov_sensors(isens) + pfso_1
502 numAss_sensors_loc(isens) = numAss_sensors_loc(isens) + 1
503 end if
504 end do
505 end if
506 end do
507
508 do familyIndex = 1, numFamily
509 call mmpi_allreduce_sumreal8scalar(tfso(familyIndex),'GRID')
510 totFSO = totFSO + tfso(familyIndex)
511 call rpn_comm_allreduce(numAss_local(familyIndex), numAss_global(familyIndex) ,1,'MPI_INTEGER','MPI_SUM','GRID',ierr)
512 end do
513
514 do isens = 1, tvs_nsensors
515 call mmpi_allreduce_sumreal8scalar(tfsotov_sensors(isens),'GRID')
516 call rpn_comm_allreduce(numAss_sensors_loc(isens), numAss_sensors_glb(isens) ,1,'MPI_INTEGER','MPI_SUM','GRID',ierr)
517 end do
518
519 if (mmpi_myid == 0) then
520
521 write(*,*) ' '
522 write(*,'(a15,f15.8)') 'Total FSO=', totFSO
523 write(*,*) ' '
524
525 do familyIndex = 1, numFamily
526 write(*,'(a4,a2,a2,f15.8,a16,i10)') 'FSO-', familyList(familyIndex), '=', tfso(familyIndex),' Count Number=', numAss_global(familyIndex)
527 end do
528 write(*,*) ' '
529
530 if (tvs_nsensors > 0) then
531 write(*,'(1x,a)') 'For TOVS decomposition by sensor:'
532 write(*,'(1x,a)') '# plt sat ins FSO'
533 do isens = 1, tvs_nsensors
534 write(*,'(i2,1x,a,1x,a,1x,i2,1x,f15.8,i10)') isens,inst_name(tvs_instruments(isens)), &
535 platform_name(tvs_platforms(isens)),tvs_satellites(isens),tfsotov_sensors(isens), numAss_sensors_glb(isens)
536 end do
537 write(*,*) ' '
538 end if
539
540 end if
541
542 end subroutine sumFSO
543
544 !--------------------------------------------------------------------------
545 ! simvar
546 !--------------------------------------------------------------------------
547 subroutine simvar(indic,nvadim,zhat,Jtotal,gradJ)
548 !
549 !:Purpose: Implement the Variational solver as described in
550 ! Courtier, 1997, Dual formulation of four-dimentional variational
551 ! assimilation, Q.J.R., pp2449-2461.
552 !
553 !:Arguments:
554 ! :indic: Value of indic
555 ! Note: 1 and 4 are reserved values for call back from m1qn3.
556 ! For direct calls use other value than 1 and 4.
557 ! =1 No action taken; =4 Both J(u) and its gradient are computed.
558 ! =2 Same as 4 (compute J and gradJ) but do not interrupt timer
559 ! of the minimizer.
560 ! =3 Compute Jo and gradJo only.
561 !
562 ! :nvadim: Dimension of the control vector in forecast error covariances space
563 !
564 ! :zhat: Control variable in forecast error covariances space
565 !
566 ! :Jtotal: Cost function of the Variational algorithm
567 !
568 ! :gradJ: Gradient of the Variational Cost funtion
569 !
570 implicit none
571
572 ! Arguments:
573 integer, intent(in) :: indic
574 integer, intent(in) :: nvadim
575 real(8), dimension(nvadim), intent(inout) :: zhat
576 real(8), intent(out) :: Jtotal
577 real(8), dimension(nvadim), intent(out) :: gradJ
578
579 ! Locals:
580 real(8) :: ahat_vhat(nvadim)
581 real(8) :: Jb, Jobs
582 type(struct_gsv) :: statevector
583 type(struct_vco), pointer :: vco_anl
584
585 call utl_tmg_stop(91)
586 call utl_tmg_stop(90)
587
588 if (indic /= 1) then ! No action taken if indic == 1
589 fso_nsim = fso_nsim + 1
590
591 if (mmpi_myid == 0) then
592 write(*,*) 'simvar: entering for simulation ',fso_nsim
593 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
594 call flush(6)
595 end if
596
597 ! note: vhat = B_t^T/2 hat(del x_t)
598 ahat_vhat(1:nvadim_mpilocal) = zhat(1:nvadim_mpilocal) + vhat(1:nvadim_mpilocal)
599
600 ! Computation of background term of cost function:
601 Jb = dot_product(zhat(1:nvadim_mpilocal),zhat(1:nvadim_mpilocal))/2.d0
602 call mmpi_allreduce_sumreal8scalar(Jb,'GRID')
603
604 vco_anl => col_getVco(columnTrlOnAnlIncLev_ptr)
605 call gsv_allocate(statevector,tim_nstepobsinc, hco_anl, vco_anl, &
606 dataKind_opt=pre_incrReal, mpi_local_opt=.true.)
607
608 call bmat_sqrtB(ahat_vhat,nvadim_mpilocal,statevector)
609
610 call s2c_tl(statevector,column_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr) ! put in column H_horiz dx
611 call utl_tmg_start(10,'--Observations')
612 call utl_tmg_start(18,'----ObsOper_TL')
613 call oop_Htl(column_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr,fso_nsim) ! Save as OBS_WORK: H_vert H_horiz dx = Hdx
614
615 call rmat_RsqrtInverseAllObs(obsSpaceData_ptr,OBS_WORK,OBS_WORK) ! Save as OBS_WORK : R**-1/2 (Hdx)
616 call utl_tmg_stop(18)
617 call utl_tmg_stop(10)
618
619 call cfn_calcJo(obsSpaceData_ptr) ! Store J-obs in OBS_JOBS : 1/2 * R**-1 (Hdx)**2
620
621 Jobs = 0.d0
622 call utl_tmg_start(90,'--Minimization')
623 call utl_tmg_start(92,'----SumCostFunction')
624 call cfn_sumJo(obsSpaceData_ptr,Jobs)
625 Jtotal = Jb + Jobs
626 if (indic == 3) then
627 Jtotal = Jobs
628 if (mmpi_myid == 0) write(*,FMT='(6X,"SIMVAR: JO = ",G23.16,6X)') Jobs
629 else
630 Jtotal = Jb + Jobs
631 if (mmpi_myid == 0) write(*,FMT='(6X,"SIMVAR: Jb = ",G23.16,6X,"JO = ",G23.16,6X,"Jt = ",G23.16)') Jb,Jobs,Jtotal
632 end if
633 call utl_tmg_stop(92)
634 call utl_tmg_stop(90)
635
636 call utl_tmg_start(10,'--Observations')
637 call rmat_RsqrtInverseAllObs(obsSpaceData_ptr,OBS_WORK,OBS_WORK) ! Modify OBS_WORK : R**-1 (Hdx)
638 call utl_tmg_stop(10)
639
640 call col_zero(column_ptr)
641
642 call utl_tmg_start(10,'--Observations')
643 call utl_tmg_start(19,'----ObsOper_AD')
644 call oop_Had(column_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr) ! Put in column : H_vert**T R**-1 (Hdx)
645 call utl_tmg_stop(19)
646 call utl_tmg_stop(10)
647
648 call s2c_ad(statevector,column_ptr,columnTrlOnAnlIncLev_ptr,obsSpaceData_ptr) ! Put in statevector H_horiz**T H_vert**T R**-1 (Hdx)
649
650 gradJ(:) = 0.d0
651 call bmat_sqrtBT(gradJ,nvadim_mpilocal,statevector)
652 call gsv_deallocate(statevector)
653
654 if (indic /= 3) then
655 gradJ(1:nvadim_mpilocal) = zhat(1:nvadim_mpilocal) + gradJ(1:nvadim_mpilocal)
656 end if
657 end if
658
659 call utl_tmg_start(90,'--Minimization')
660 call utl_tmg_start(91,'----QuasiNewton')
661
662 if (mmpi_myid == 0) write(*,*) 'simvar: Finished'
663
664 end subroutine simvar
665
666 !--------------------------------------------------------------------------
667 ! prscal
668 !--------------------------------------------------------------------------
669 subroutine prscal(kdim,px,py,ddsc)
670 !
671 !:Purpose: evaluation of the inner product in canonical space used in the minimization
672 !
673 implicit none
674
675 ! Arguments:
676 integer, intent(in) :: kdim ! dimension of the vectors
677 real(8), intent(in) :: px(kdim) ! vector components for which <px,py> is being calculated
678 real(8), intent(in) :: py(kdim) ! vector components for which <px,py> is being calculated
679 real(8), intent(out) :: ddsc ! result of the inner product
680
681 ! Locals:
682 integer :: cvIndex
683
684 ddsc = 0.d0
685
686 do cvIndex=1,nvadim_mpilocal
687 ddsc = ddsc + px(cvIndex)*py(cvIndex)
688 end do
689
690 call mmpi_allreduce_sumreal8scalar(ddsc,'GRID')
691
692 end subroutine prscal
693
694 !--------------------------------------------------------------------------
695 ! dcanab
696 !--------------------------------------------------------------------------
697 subroutine dcanab(kdim,py,px)
698 !
699 !:Purpose: Change of variable associated with the canonical inner product
700 ! to compute PX = L^-1 * Py with L related to the inner product
701 ! <PX,PY> = PX^t L^t L PY
702 ! (see the modulopt documentation aboutn DTCAB)
703 ! Double precision version based on single precision CTCAB.
704 ! Refered to as dummy argument DTCAB by N1QN3 minimization
705 ! package.
706 !
707 implicit none
708
709 ! Arguments:
710 integer, intent(in) :: kdim ! dimension of the vectors
711 real(8), intent(inout) :: px(kdim)
712 real(8), intent(in) :: py(kdim)
713
714 ! Locals:
715 integer :: cvIndex
716
717 do cvIndex = 1, kdim
718 px(cvIndex) = py(cvIndex)
719 end do
720
721 end subroutine dcanab
722
723 !--------------------------------------------------------------------------
724 ! dcanonb
725 !--------------------------------------------------------------------------
726 subroutine dcanonb(kdim,px,py)
727 !
728 !:Purpose: Change of variable associated with the canonical inner product
729 ! to compute PY = L * PX with L related to the inner product
730 ! <PX,PY> = PX^t L^t L PY
731 ! (see the modulopt documentation about DTONB)
732 ! Double precision version based on single precision CTCAB.
733 ! Refered to as dummy argument DTCAB by N1QN3 minimization
734 ! package.
735 !
736 implicit none
737
738 ! Arguments:
739 integer, intent(in) :: kdim ! dimension of the vectors
740 real(8), intent(in) :: px(kdim)
741 real(8), intent(inout) :: py(kdim)
742
743 ! Locals:
744 integer :: cvIndex
745
746 do cvIndex = 1, kdim
747 py(cvIndex) = px(cvIndex)
748 end do
749
750 end subroutine dcanonb
751
752 !--------------------------------------------------------------------------
753 ! multEnergyNorm
754 !--------------------------------------------------------------------------
755 subroutine multEnergyNorm(statevector_inout, statevector_ref, &
756 latMin, latMax, lonMin, lonMax, &
757 uvNorm,ttNorm,p0Norm,huNorm,tgNorm, straNorm)
758 !
759 !:Purpose: Computes energy norms
760 ! For some positive definite symmetric matrix defining the energy,
761 ! total energy = x^T * C * x
762 ! statevector_inout = C * statevector_inout
763 ! (Buehner, Du and Bedard, 2018)
764 !
765 implicit none
766
767 ! Arguments:
768 type(struct_gsv), intent(inout) :: statevector_inout
769 type(struct_gsv), intent(in) :: statevector_ref
770 real(8), intent(in) :: latMin
771 real(8), intent(in) :: latMax
772 real(8), intent(in) :: lonMin
773 real(8), intent(in) :: lonMax
774 logical, intent(in) :: uvNorm
775 logical, intent(in) :: ttNorm
776 logical, intent(in) :: p0Norm
777 logical, intent(in) :: huNorm
778 logical, intent(in) :: tgNorm
779 logical, intent(in) :: straNorm
780
781 ! Locals:
782 integer :: stepIndex, lonIndex, levIndex, latIndex, lonIndex2, latIndex2, nLev_M, nLev_T
783 real(8) :: scaleFactor, scaleFactorConst, scaleFactorLat, scaleFactorLon, scaleFactorLev
784 real(8) :: pfac, tfac, qfac
785 real(8) :: sumScale , sumeu, sumev, sumep, sumet, sumeq
786 real(8), pointer :: field_UU(:,:,:,:), field_VV(:,:,:,:), field_T(:,:,:,:), field_LQ(:,:,:,:)
787 real(8), pointer :: field_Psfc(:,:,:,:), field_TG(:,:,:,:),Psfc_ptr(:,:,:)
788 real(8), pointer :: Press_T(:,:,:)
789 real(8), pointer :: Press_M(:,:,:)
790 real(8), allocatable :: Psfc_ref(:,:)
791 real(8), parameter :: T_r = 280.0D0
792 real(8), parameter :: Psfc_r = 100000.0D0 ! unit Pa
793 real(8), parameter :: sigma = 0.3 ! weight factor for humidity
794
795 if (mmpi_myid == 0) write(*,*) 'multEnergyNorm: START'
796 nullify(Press_T,Press_M)
797
798 ! the factors for TT, HU and Ps (for wind is 1)
799 tfac = mpc_cp_dry_air_r8/T_r ! temperature factor (c_p/T_r)
800 qfac = sigma*mpc_heat_condens_water_r8**2/(mpc_cp_dry_air_r8*T_r) ! humidity factor ( (l_p*l_p)/(c_p*T_r) )
801 pfac = mpc_rgas_dry_air_r8*T_r/(Psfc_r**2) ! surface pressure factor (R*T_r/Psfc_r^2)
802
803 if (.not. gsv_isAllocated(statevector_inout)) then
804 call utl_abort('multEnergyNorm: gridStateVector_inout not yet allocated')
805 end if
806
807 nLev_M = gsv_getNumLev(statevector_inout,'MM')
808 nLev_T = gsv_getNumLev(statevector_inout,'TH')
809
810 ! compute 3D log pressure fields
811 call gsv_getField(statevector_ref,Psfc_ptr,'P0')
812 allocate(Psfc_ref(statevector_inout%lonPerPEmax,statevector_inout%latPerPEmax))
813 Psfc_ref(:,:) = &
814 Psfc_ptr(statevector_inout%myLonBeg:statevector_inout%myLonEnd, &
815 statevector_inout%myLatBeg:statevector_inout%myLatEnd, 1)
816 call czp_fetch3DLevels(statevector_inout%vco, Psfc_ref, &
817 fldM_opt=Press_M, fldT_opt=Press_T)
818 ! dlat * dlon
819 scaleFactorConst = statevector_inout%hco%dlat*statevector_inout%hco%dlon
820
821 ! for wind components if to include in Norm calculation
822 call gsv_getField(statevector_inout,field_UU,'UU')
823 call gsv_getField(statevector_inout,field_VV,'VV')
824
825 sumeu = 0.0D0
826 sumev = 0.0D0
827 sumScale = 0.0D0
828 if (uvNorm) then
829 do levIndex = 1, nLev_M
830 do stepIndex = 1, statevector_inout%numStep
831 do latIndex = statevector_inout%myLatBeg, statevector_inout%myLatEnd
832 latIndex2 = latIndex - statevector_inout%myLatBeg + 1
833 ! IF lat is out of the domain where we want to compute the NRJ norm, we put scaleFactorLat = 0.
834 if (statevector_inout%hco%lat(latIndex) >= latMin .and. statevector_inout%hco%lat(latIndex) <= latMax) then
835 scaleFactorLat = cos(statevector_inout%hco%lat(latIndex))
836 else
837 scaleFactorLat = 0.0D0
838 end if
839 do lonIndex = statevector_inout%myLonBeg, statevector_inout%myLonEnd
840 lonIndex2 = lonIndex - statevector_inout%myLonBeg + 1
841 ! Similarly, if lon is out of the domain where we want to compute the NRJ norm, we put scaleFactorLon = 0.
842 if (statevector_inout%hco%lon(lonIndex) >= lonMin .and. statevector_inout%hco%lon(lonIndex) <= lonMax) then
843 scaleFactorLon = 1.0D0
844 else
845 scaleFactorLon = 0.0D0
846 end if
847 ! do all thermo levels for which there is a momentum level above and below
848 if (straNorm) then !for strato Norm
849 if (levIndex == nLev_M .or. levIndex == 1 .or. &
850 Press_T(lonIndex2, latIndex2, levIndex) < 100.0D0 .or. &
851 Press_T(lonIndex2, latIndex2, levIndex) > 10000.0D0) then
852 scaleFactorLev = 0.0D0
853 else
854 scaleFactorLev = Press_T(lonIndex2, latIndex2, levIndex+1) - Press_T(lonIndex2, latIndex2, levIndex)
855 end if
856
857 else
858 if ( levIndex == nLev_M) then
859 scaleFactorLev = Press_M(lonIndex2, latIndex2, nLev_M)-Press_T(lonIndex2, latIndex2, nLev_T-1)
860 else if ( Press_T(lonIndex2, latIndex2, levIndex) < 10000.0D0) then
861 scaleFactorLev = 0.0D0
862 else
863 scaleFactorLev = Press_T(lonIndex2, latIndex2, levIndex+1) - Press_T(lonIndex2, latIndex2, levIndex)
864 end if
865 end if
866
867 scaleFactor = scaleFactorConst * scaleFactorLat* scaleFactorLon * scaleFactorLev
868 sumScale = sumScale + scaleFactor
869
870 sumeu = sumeu + &
871 0.5 * field_UU(lonIndex,latIndex,levIndex,stepIndex) * field_UU(lonIndex,latIndex,levIndex,stepIndex) * scaleFactor
872 sumev = sumev + &
873 0.5 * field_VV(lonIndex,latIndex,levIndex,stepIndex) * field_VV(lonIndex,latIndex,levIndex,stepIndex) * scaleFactor
874
875 field_UU(lonIndex,latIndex,levIndex,stepIndex) = &
876 field_UU(lonIndex,latIndex,levIndex,stepIndex) * 0.5 * scaleFactor
877 field_VV(lonIndex,latIndex,levIndex,stepIndex) = &
878 field_VV(lonIndex,latIndex,levIndex,stepIndex) * 0.5 * scaleFactor
879 end do !lonIndex
880 end do !latIndex
881 end do ! stepIndex
882 end do ! levIndex
883 call mmpi_allreduce_sumreal8scalar(sumeu,'grid')
884 call mmpi_allreduce_sumreal8scalar(sumev,'grid')
885 call mmpi_allreduce_sumreal8scalar(sumScale,'grid')
886 sumeu = sumeu/sumScale
887 sumev = sumev/sumScale
888
889 field_UU(:,:,:,:) = field_UU(:,:,:,:)/sumScale
890 field_VV(:,:,:,:) = field_VV(:,:,:,:)/sumScale
891 else
892 field_UU(:,:,:,:) = field_UU(:,:,:,:)*0.0D0
893 field_VV(:,:,:,:) = field_VV(:,:,:,:)*0.0D0
894 end if ! if uvNorm
895
896 if (mmpi_myid == 0) write(*,*) 'energy for UU=', sumeu
897 if (mmpi_myid == 0) write(*,*) 'energy for VV=', sumev
898
899 ! for Temperature
900 call gsv_getField(statevector_inout,field_T,'TT')
901 sumScale = 0.0D0
902 sumet = 0.0D0
903 if (ttNorm) then
904 do levIndex = 1, nLev_T
905 do stepIndex = 1, statevector_inout%numStep
906 do latIndex = statevector_inout%myLatBeg, statevector_inout%myLatEnd
907 latIndex2 = latIndex - statevector_inout%myLatBeg + 1
908 ! IF lat is out of the domain where we want to compute the NRJ norm, we put scaleFactorLat = 0.
909 if (statevector_inout%hco%lat(latIndex) >= latMin .and. statevector_inout%hco%lat(latIndex) <= latMax) then
910 scaleFactorLat = cos(statevector_inout%hco%lat(latIndex))
911 else
912 scaleFactorLat = 0.0D0
913 end if
914 do lonIndex = statevector_inout%myLonBeg, statevector_inout%myLonEnd
915 lonIndex2 = lonIndex - statevector_inout%myLonBeg + 1
916 ! Similarly, if lon is out of the domain where we want to compute the NRJ norm, we put scaleFactorLon = 0.
917 if (statevector_inout%hco%lon(lonIndex) >= lonMin .and. statevector_inout%hco%lon(lonIndex) <= lonMax) then
918 scaleFactorLon = 1.0D0
919 else
920 scaleFactorLon = 0.0D0
921 end if
922 ! do all thermo levels for which there is a momentum level above and below
923 if (straNorm) then !for strato norm
924 if (levIndex == nLev_T .or. levIndex == 1) then
925 scaleFactorLev = 0.0D0
926 else if (Press_M(lonIndex2, latIndex2, levIndex-1) < 100.0D0 .or. &
927 Press_M(lonIndex2, latIndex2, levIndex-1) > 10000.0D0 ) then
928 scaleFactorLev = 0.0D0
929 else
930 scaleFactorLev = Press_M(lonIndex2, latIndex2, levIndex ) - Press_M(lonIndex2, latIndex2, levIndex-1)
931 end if
932
933 else
934 if (levIndex == nLev_T) then !surface
935 scaleFactorLev = Press_T(lonIndex2, latIndex2, nLev_T)-Press_T(lonIndex2, latIndex2, nLev_T-1)
936 else if (levIndex == 1) then ! top
937 scaleFactorLev = 0.0D0
938 else if ( Press_M(lonIndex2, latIndex2, levIndex-1) < 10000.0D0) then
939 scaleFactorLev = 0.0D0
940 else
941 scaleFactorLev = Press_M(lonIndex2, latIndex2, levIndex ) - Press_M(lonIndex2, latIndex2, levIndex-1)
942 end if
943 end if
944
945 scaleFactor = scaleFactorConst * scaleFactorLat * scaleFactorLon * scaleFactorLev
946 sumet = sumet + &
947 0.5 * tfac * field_T(lonIndex,latIndex,levIndex,stepIndex) * field_T(lonIndex,latIndex,levIndex,stepIndex) * scaleFactor
948 sumScale = sumScale + scaleFactor
949 field_T(lonIndex,latIndex,levIndex,stepIndex) = &
950 field_T(lonIndex,latIndex,levIndex,stepIndex) * 0.5 * tfac * scaleFactor
951 end do
952 end do
953 end do ! stepIndex
954 end do ! levIndex
955 call mmpi_allreduce_sumreal8scalar(sumet,'grid')
956 call mmpi_allreduce_sumreal8scalar(sumScale,'grid')
957 sumet = sumet/sumScale
958 field_T(:,:,:,:) = field_T(:,:,:,:)/sumScale
959 else
960 field_T(:,:,:,:) = field_T(:,:,:,:)*0.0D0
961 end if ! if ttNorm
962
963 if (mmpi_myid == 0) write(*,*) 'energy for TT=', sumet
964
965 ! humidity (set to zero, for now)
966 call gsv_getField(statevector_inout,field_LQ,'HU')
967 sumScale = 0.0D0
968 sumeq = 0.0D0
969 if (huNorm) then
970 do levIndex = 1, nLev_T
971 do stepIndex = 1, statevector_inout%numStep
972 do latIndex = statevector_inout%myLatBeg, statevector_inout%myLatEnd
973 latIndex2 = latIndex - statevector_inout%myLatBeg + 1
974 ! IF lat is out of the domain where we want to compute the NRJ norm, we put scaleFactorLat = 0.
975 if (statevector_inout%hco%lat(latIndex) >= latMin .and. statevector_inout%hco%lat(latIndex) <= latMax) then
976 scaleFactorLat = cos(statevector_inout%hco%lat(latIndex))
977 else
978 scaleFactorLat = 0.0D0
979 end if
980 do lonIndex = statevector_inout%myLonBeg, statevector_inout%myLonEnd
981 lonIndex2 = lonIndex - statevector_inout%myLonBeg + 1
982 ! Similarly, if lon is out of the domain where we want to compute the NRJ norm, we put scaleFactorLon = 0.
983 if (statevector_inout%hco%lon(lonIndex) >= lonMin .and. statevector_inout%hco%lon(lonIndex) <= lonMax) then
984 scaleFactorLon = 1.0D0
985 else
986 scaleFactorLon = 0.0D0
987 end if
988 ! do all thermo levels for which there is a momentum level above and below
989 if (straNorm) then !for strato norm
990 if (levIndex == nLev_T .or. levIndex == 1) then
991 scaleFactorLev = 0.0D0
992 else if(Press_M(lonIndex2, latIndex2, levIndex-1) < 100.0D0 .or. &
993 Press_M(lonIndex2, latIndex2, levIndex-1) > 10000.0D0 ) then
994 scaleFactorLev = 0.0D0
995 else
996 scaleFactorLev = Press_M(lonIndex2, latIndex2, levIndex ) - Press_M(lonIndex2, latIndex2, levIndex-1)
997 end if
998
999 else
1000 if ( levIndex == nLev_T) then !surface
1001 scaleFactorLev = Press_T(lonIndex2, latIndex2, nLev_T) - Press_T(lonIndex2, latIndex2, nLev_T-1)
1002 else if (levIndex == 1) then ! top
1003 scaleFactorLev = 0.0D0
1004 else if ( Press_M(lonIndex2, latIndex2, levIndex-1) < 10000.0D0) then
1005 scaleFactorLev = 0.0D0
1006 else
1007 scaleFactorLev = Press_M(lonIndex2, latIndex2, levIndex ) - Press_M(lonIndex2, latIndex2, levIndex-1)
1008 end if
1009 end if
1010
1011 scaleFactor = scaleFactorConst * scaleFactorLat * scaleFactorLon * scaleFactorLev
1012 sumScale = sumScale + scaleFactor
1013
1014 sumeq = sumeq + 0.5 * qfac * &
1015 field_LQ(lonIndex,latIndex,levIndex,stepIndex) * field_LQ(lonIndex,latIndex,levIndex,stepIndex) * scaleFactor
1016
1017 field_LQ(lonIndex,latIndex,levIndex,stepIndex) = &
1018 field_LQ(lonIndex,latIndex,levIndex,stepIndex) * 0.5 * scaleFactor * qfac
1019
1020 end do
1021 end do
1022 end do ! stepIndex
1023 end do ! latIndex
1024 call mmpi_allreduce_sumreal8scalar(sumScale,'grid')
1025 call mmpi_allreduce_sumreal8scalar(sumeq,'grid')
1026 sumeq = sumeq/sumScale
1027 field_LQ(:,:,:,:) = field_LQ(:,:,:,:)/sumScale
1028 else
1029 field_LQ(:,:,:,:) = field_LQ(:,:,:,:)*0.0D0
1030 end if ! if huNorm
1031
1032 if (mmpi_myid == 0) write(*,*) 'energy for HU=', sumeq
1033
1034 ! surface pressure
1035 call gsv_getField(statevector_inout,field_Psfc,'P0')
1036 sumScale = 0.0D0
1037 sumep = 0.0
1038 if (p0Norm .and. .not.straNorm) then
1039 do stepIndex = 1, statevector_inout%numStep
1040 do latIndex = statevector_inout%myLatBeg, statevector_inout%myLatEnd
1041 ! IF lat is out of the domain where we want to compute the NRJ norm, we put scaleFactorLat = 0.
1042 if (statevector_inout%hco%lat(latIndex) >= latMin .and. statevector_inout%hco%lat(latIndex) <= latMax) then
1043 scaleFactorLat = cos(statevector_inout%hco%lat(latIndex))
1044 else
1045 scaleFactorLat = 0.0D0
1046 end if
1047 do lonIndex = statevector_inout%myLonBeg, statevector_inout%myLonEnd
1048 ! Similarly, if lon is out of the domain where we want to compute the NRJ norm, we put scaleFactorLon = 0.
1049 if (statevector_inout%hco%lon(lonIndex) >= lonMin .and. statevector_inout%hco%lon(lonIndex) <= lonMax) then
1050 scaleFactorLon = 1.0D0
1051 else
1052 scaleFactorLon = 0.0D0
1053 end if
1054 scaleFactor = scaleFactorConst * scaleFactorLat * scaleFactorLon
1055 sumScale = sumScale + scaleFactor
1056 sumep = sumep + 0.5 * pfac * &
1057 field_Psfc(lonIndex,latIndex,1,stepIndex) * field_Psfc(lonIndex,latIndex,1,stepIndex) * scaleFactor
1058 field_Psfc(lonIndex,latIndex,1,stepIndex) = &
1059 field_Psfc(lonIndex,latIndex,1,stepIndex) * 0.5 * scaleFactor * pfac
1060 end do
1061 end do ! latIndex
1062 end do ! stepIndex
1063
1064 call mmpi_allreduce_sumreal8scalar(sumep,'grid')
1065 call mmpi_allreduce_sumreal8scalar(sumScale,'grid')
1066 sumep = sumep/sumScale
1067 field_Psfc(:,:,:,:) = field_Psfc(:,:,:,:)/sumScale
1068 else
1069 field_Psfc(:,:,:,:) = field_Psfc(:,:,:,:)*0.0D0
1070 end if ! if p0Norm
1071
1072 if (mmpi_myid == 0) write(*,*) 'energy for Ps=', sumep
1073
1074 ! skin temperature (set to zero for now)
1075 call gsv_getField(statevector_inout,field_TG,'TG')
1076 sumScale = 0.0D0
1077 if (tgNorm .and. .not.straNorm) then
1078 do stepIndex = 1, statevector_inout%numStep
1079 do latIndex = statevector_inout%myLatBeg, statevector_inout%myLatEnd
1080 ! IF lat is out of the domain where we want to compute the NRJ norm, we put scaleFactorLat = 0.
1081 if (statevector_inout%hco%lat(latIndex) >= latMin .and. statevector_inout%hco%lat(latIndex) <= latMax) then
1082 scaleFactorLat = cos(statevector_inout%hco%lat(latIndex))
1083 else
1084 scaleFactorLat = 0.0D0
1085 end if
1086 do lonIndex = statevector_inout%myLonBeg, statevector_inout%myLonEnd
1087 ! Similarly, if lon is out of the domain where we want to compute the NRJ norm, we put scaleFactorLon = 0.
1088 if (statevector_inout%hco%lon(lonIndex) >= lonMin .and. statevector_inout%hco%lon(lonIndex) <= lonMax) then
1089 scaleFactorLon = 1.0D0
1090 else
1091 scaleFactorLon = 0.0D0
1092 end if
1093 scaleFactor = scaleFactorConst * scaleFactorLat * scaleFactorLon
1094 sumScale = sumScale + scaleFactor
1095 field_TG(lonIndex,latIndex,1,stepIndex) = &
1096 field_TG(lonIndex,latIndex,1,stepIndex) * 0.5 * scaleFactor * 0.0
1097 end do
1098 end do ! latIndex
1099 end do ! stepIndex
1100 call mmpi_allreduce_sumreal8scalar(sumScale,'grid')
1101 field_TG(:,:,:,:) = field_TG(:,:,:,:)/sumScale
1102 else
1103 field_TG(:,:,:,:) = field_TG(:,:,:,:)*0.0D0
1104 end if ! if tgNorm
1105
1106 if (mmpi_myid == 0) write(*,*) 'energy for total=', sumeu + sumev + sumet + sumep + sumeq
1107 deallocate(Press_T,Press_M)
1108 deallocate(Psfc_ref)
1109
1110 if (mmpi_myid == 0) write(*,*) 'multEnergyNorm: END'
1111
1112 end subroutine multEnergyNorm
1113
1114end module fsoi_mod
