1module oceanBackground_mod
2 ! MODULE oceanBackground_mod (prefix='obgd' category='1. High-level functionality')
3 !
4 !:Purpose: storage for ocean background related subroutines
5 !
6 use codtyp_mod
7 use gridStateVector_mod
8 use gridStateVectorFileIO_mod
9 use horizontalCoord_mod
10 use verticalCoord_mod
11 use timeCoord_mod
12
13 implicit none
14
15 save
16 private
17
18 ! Public functions/subroutines
19 public :: obgd_computeSSTrial
20
21 ! External functions
22 integer, external :: fnom, fclos
23
24 contains
25
26 !----------------------------------------------------------------------------------------
27 ! obgd_computeSSTrial
28 !----------------------------------------------------------------------------------------
29 subroutine obgd_computeSSTrial(hco, vco, trialDateStamp, analysisDateStamp, &
30 nmonthsClim, datestampClim, alphaClim, etiket)
31 !
32 !: Purpose: to compute SST background
33 ! xb(t) = (xa(t-1) - xclim(t-1))*alpha + xclim(t)
34 implicit none
35
36 ! Arguments:
37 type(struct_hco), pointer, intent(in) :: hco ! horizontal grid
38 type(struct_vco), pointer, intent(in) :: vco ! vertical grid
39 integer , intent(in) :: trialDateStamp ! trial datestamp
40 integer , intent(in) :: analysisDateStamp ! datestamp for last analysis
41 integer , intent(in) :: nmonthsClim ! number of climatological fields (= 12)
42 integer , intent(in) :: datestampClim(:) ! datestamps of input climatology fields
43 real(4) , intent(in) :: alphaClim ! scalling factor to relax towards climatology
44 character(len=10) , intent(in) :: etiket ! etiket from namelist and for trial
45
46 ! Locals:
47 type(struct_gsv) :: stateVector
48 real(4), pointer :: stateVector_ptr(:, :, :)
49 integer :: lonIndex, latIndex
50 real(8) :: climatology_m1(hco % ni, hco % nj), climatology(hco % ni, hco % nj)
51
52 write(*,*) 'obgd_computeSSTrial: starting...'
53
54 ! get SST analysis
55 call gsv_allocate(stateVector, 1, hco, vco, dataKind_opt = 4, datestamp_opt = analysisDateStamp, &
56 mpi_local_opt = .true., varNames_opt = (/'TM'/))
57 call gio_readFromFile(stateVector, './analysis', ' ','A', &
58 unitConversion_opt=.false., containsFullField_opt=.true.)
59 call gsv_getField(stateVector, stateVector_ptr)
60
61 call obgd_getClimatology(analysisDateStamp, hco, vco, nmonthsClim, datestampClim, climatology_m1)
62 call obgd_getClimatology(trialDateStamp , hco, vco, nmonthsClim, datestampClim, climatology)
63
64 ! compute background state
65 ! xb(t) = (xa(t-1) - xclim(t-1))*alpha + xclim(t)
66 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
67 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
68 stateVector_ptr(lonIndex, latIndex, 1) = (stateVector_ptr(lonIndex, latIndex, 1) - &
69 climatology_m1(lonIndex, latIndex)) * &
70 alphaClim + climatology(lonIndex, latIndex)
71 end do
72 end do
73
74 ! modify dateStamp (from analysis) with trial dateStamp
75 call gsv_modifyDate(stateVector, trialDateStamp, modifyDateOrigin_opt = .true.)
76
77 ! save trial field
78 call gio_writeToFile(stateVector, './trial', etiket, typvar_opt = 'P@')
79
80 call gsv_deallocate(stateVector)
81
82 end subroutine obgd_computeSSTrial
83
84 !----------------------------------------------------------------------------------------
85 ! obgd_getClimatology
86 !----------------------------------------------------------------------------------------
87 subroutine obgd_getClimatology(dateStamp, hco, vco, nmonthsClim, datestampClim, output)
88 !
89 !: Purpose: 1) to read SST climatological fields from a std file
90 ! 2) to interpolate the field in time using the current day (t) in current month (m)
91 ! SST(t) = SST_clim(m) + (t-1)/(ndays-1) * (SST_clim(m+1) - SST_clim(m)),
92 ! where ndays is a number of days in current month
93 !
94 implicit none
95
96 ! Arguments:
97 integer , intent(in) :: dateStamp ! date stamp for the current day
98 type(struct_hco), pointer, intent(in) :: hco ! horizontal grid
99 type(struct_vco), pointer, intent(in) :: vco ! vertical grid
100 integer , intent(in) :: nmonthsClim ! number of records in the climatology file
101 integer , intent(in) :: datestampClim(:) ! datestamps in the climatology file
102 real(8) , intent(inout) :: output(:,:) ! interpolated SST field from climatology
103
104 ! Locals:
105 integer :: hour, day, month, yyyy, ndays, nextMonth
106 type(struct_gsv) :: stateVector, stateVector_nextMonth
107 real(4), pointer :: clim_ptr(:, :, :), clim_nextMonth_ptr(:, :, :)
108 integer :: lonIndex, latIndex
109
110 call tim_dateStampToYYYYMMDDHH(dateStamp, hour, day, month, ndays, yyyy)
111 write(*,'(a,3i5,a,i12,a)') 'obgd_getClimatology: interpolating climatology for day/month/year (datestamp): ', &
112 day, month, yyyy, '(', datestamp, ')'
113
114 if (month == nmonthsClim) then
115 nextMonth = 1
116 else
117 nextMonth = month + 1
118 end if
119
120 ! get climatology, current month
121 write(*,*) 'obgd_getClimatology: reading climatology, month: ', month, ', datestamp: ', datestampClim(month)
122 call gsv_allocate(stateVector, 1, hco, vco, dataKind_opt = 4, &
123 datestamp_opt = datestampClim(month), mpi_local_opt = .true., &
124 varNames_opt = (/'TM'/), hInterpolateDegree_opt ='LINEAR')
125 call gio_readFromFile(stateVector, './climatology', ' ',' ', &
126 unitConversion_opt=.false., containsFullField_opt=.true.)
127 call gsv_getField(stateVector, clim_ptr)
128
129 ! get climatology, next month
130 write(*,*) 'obgd_getClimatology: reading climatology, month: ', nextMonth, ', datestamp: ', datestampClim(nextMonth)
131 call gsv_allocate(stateVector_nextMonth, 1, hco, vco, dataKind_opt = 4, &
132 datestamp_opt = datestampClim(nextMonth), mpi_local_opt = .true., &
133 varNames_opt = (/'TM'/), hInterpolateDegree_opt ='LINEAR')
134 call gio_readFromFile(stateVector_nextMonth, './climatology', ' ',' ', &
135 unitConversion_opt=.false., containsFullField_opt=.true.)
136 call gsv_getField(stateVector_nextMonth, clim_nextMonth_ptr)
137
138 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
139 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
140 output(lonIndex, latIndex) = clim_ptr(lonIndex, latIndex, 1) + real(day - 1) / real(ndays - 1) * &
141 (clim_nextMonth_ptr(lonIndex, latIndex, 1) - clim_ptr(lonIndex, latIndex, 1))
142 end do
143 end do
144
145 call gsv_deallocate(stateVector)
146 call gsv_deallocate(stateVector_nextMonth)
147
148 end subroutine obgd_getClimatology
149
150end module oceanBackground_mod
