/[MITgcm]/MITgcm/pkg/ecco/cost_ctdt.F
ViewVC logotype

Contents of /MITgcm/pkg/ecco/cost_ctdt.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.12 - (show annotations) (download)
Sat Oct 18 18:15:44 2014 UTC (9 years, 7 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65p, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65g
Changes since 1.11: +4 -32 lines
- add CPP brackets around includes, to omit
  them altogether when they are not used.

1 C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ctdt.F,v 1.11 2012/08/10 19:45:25 jmc Exp $
2 C $Name: $
3
4 #include "ECCO_OPTIONS.h"
5
6
7 subroutine cost_CTDT(
8 I myiter,
9 I mytime,
10 I mythid
11 & )
12
13 c ==================================================================
14 c SUBROUTINE cost_CTDT
15 c ==================================================================
16 c
17 c o Evaluate cost function contribution of CTD temperature data.
18 c
19 c started: Elisabeth Remy eremy@ucsd.edu 30-Aug-2000
20 c
21 c
22 c ==================================================================
23 c SUBROUTINE cost_CTDT
24 c ==================================================================
25
26 implicit none
27
28 c == global variables ==
29
30 #ifdef ALLOW_CTDT_COST_CONTRIBUTION
31 #include "EEPARAMS.h"
32 #include "SIZE.h"
33 #include "GRID.h"
34 #include "DYNVARS.h"
35
36 #include "cal.h"
37 #include "ecco_cost.h"
38 #include "CTRL_SIZE.h"
39 #include "ctrl.h"
40 #include "ctrl_dummy.h"
41 #include "optim.h"
42 #endif
43
44 c == routine arguments ==
45
46 integer myiter
47 _RL mytime
48 integer mythid
49
50 c == local variables ==
51
52 #ifdef ALLOW_CTDT_COST_CONTRIBUTION
53 integer bi,bj
54 integer i,j,k
55 integer itlo,ithi
56 integer jtlo,jthi
57 integer jmin,jmax
58 integer imin,imax
59 integer nrec
60 integer irec
61 integer ilu
62
63 _RL fctile_ctdt
64 _RL fcthread_ctdt
65 _RL www (1-olx:snx+olx,1-oly:sny+oly)
66 _RL wtmp (1-olx:snx+olx,1-oly:sny+oly)
67 _RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
68 _RL tmpbar (1-olx:snx+olx,1-oly:sny+oly)
69 _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
70 _RL spval
71 _RL spmax
72
73 character*(80) fnametheta
74
75 logical doglobalread
76 logical ladinit
77
78 character*(MAX_LEN_MBUF) msgbuf
79
80 cnew(
81 integer il
82 integer mody, modm
83 integer iyear, imonth
84 character*(80) fnametmp
85 logical exst
86 cnew)
87
88 c == external functions ==
89
90 integer ilnblnk
91 external ilnblnk
92
93 c == end of interface ==
94
95 jtlo = mybylo(mythid)
96 jthi = mybyhi(mythid)
97 itlo = mybxlo(mythid)
98 ithi = mybxhi(mythid)
99 jmin = 1
100 jmax = sny
101 imin = 1
102 imax = snx
103
104 spval = -1.8
105 spmax = 40.
106
107 c-- Read state record from global file.
108 doglobalread = .false.
109 ladinit = .false.
110
111 if (optimcycle .ge. 0) then
112 ilu=ilnblnk( tbarfile )
113 write(fnametheta(1:80),'(2a,i10.10)')
114 & tbarfile(1:ilu),'.',optimcycle
115 endif
116
117 fcthread_ctdt = 0. _d 0
118
119 cnew(
120 mody = modelstartdate(1)/10000
121 modm = modelstartdate(1)/100 - mody*100
122 cnew)
123
124 c-- Loop over records.
125 do irec = 1,nmonsrec
126
127 c-- Read time averages and the monthly mean data.
128 call active_read_xyz( fnametheta, tbar, irec,
129 & doglobalread, ladinit,
130 & optimcycle, mythid, xx_tbar_mean_dummy )
131
132 cnew(
133 iyear = mody + INT((modm-1+irec-1)/12)
134 imonth = 1 + MOD(modm-1+irec-1,12)
135 il=ilnblnk(ctdtfile)
136 write(fnametmp(1:80),'(2a,i4)')
137 & ctdtfile(1:il), '_', iyear
138 inquire( file=fnametmp, exist=exst )
139 if (.NOT. exst) then
140 write(fnametmp(1:80),'(a)') ctdtfile(1:il)
141 imonth = irec
142 endif
143
144 call mdsreadfield( fnametmp, cost_iprec, 'RL', nr, ctdtobs,
145 & imonth, mythid)
146 cnew)
147
148 c-- Loop over this thread tiles.
149 do bj = jtlo,jthi
150 do bi = itlo,ithi
151 c-- Loop over the model layers
152
153 fctile_ctdt = 0. _d 0
154
155 do k = 1,nr
156
157 c-- Determine the weights to be used.
158 do j = jmin,jmax
159 do i = imin,imax
160 cmask(i,j) = 1. _d 0
161 if (ctdtobs(i,j,k,bi,bj) .lt. spval .or.
162 & ctdtobs(i,j,k,bi,bj) .gt. spmax .or.
163 & ctdtobs(i,j,k,bi,bj) .eq. 0. ) then
164 cmask(i,j) = 0. _d 0
165 endif
166
167 c set cmask=0 in areas shallower than 1000m
168
169 if ( _hFacC(i,j,k,bi,bj) .ne. 0. ) then
170
171 www(i,j) = cosphi(i,j,bi,bj)*cmask(i,j)
172 tmpobs(i,j) = ctdtobs(i,j,k,bi,bj)
173 tmpbar(i,j) = tbar(i,j,k,bi,bj)
174 wtmp(i,j) = wtheta2(i,j,k,bi,bj)
175
176 c-- The array ctdtobs contains CTD temperature.
177 fctile_ctdt = fctile_ctdt +
178 & (wtmp(i,j)*www(i,j))*
179 & (tmpbar(i,j)-tmpobs(i,j))*
180 & (tmpbar(i,j)-tmpobs(i,j))
181 if ( wtmp(i,j)*www(i,j) .ne. 0. )
182 & num_ctdt(bi,bj) = num_ctdt(bi,bj) + 1. _d 0
183 endif
184 enddo
185 enddo
186 enddo
187 c-- End of loop over layers.
188
189 fcthread_ctdt = fcthread_ctdt + fctile_ctdt
190 objf_ctdt(bi,bj) = objf_ctdt(bi,bj) + fctile_ctdt
191
192 #ifdef ECCO_VERBOSE
193 write(msgbuf,'(a)') ' '
194 call print_message( msgbuf, standardmessageunit,
195 & SQUEEZE_RIGHT , mythid)
196 write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
197 & ' COST_CTDT: irec,bi,bj = ',irec,bi,bj
198 call print_message( msgbuf, standardmessageunit,
199 & SQUEEZE_RIGHT , mythid)
200 write(msgbuf,'(a,d22.15)')
201 & ' COST_CTDT: cost function = ', fctile_ctdt
202 call print_message( msgbuf, standardmessageunit,
203 & SQUEEZE_RIGHT , mythid)
204 write(msgbuf,'(a)') ' '
205 call print_message( msgbuf, standardmessageunit,
206 & SQUEEZE_RIGHT , mythid)
207 #endif
208
209 enddo
210 enddo
211
212 #ifdef ECCO_VERBOSE
213 c-- Print cost function for all tiles.
214 _GLOBAL_SUM_RL( fcthread_ctdt , myThid )
215 write(msgbuf,'(a)') ' '
216 call print_message( msgbuf, standardmessageunit,
217 & SQUEEZE_RIGHT , mythid)
218 write(msgbuf,'(a,i8.8)')
219 & ' cost_CTDT: irec = ',irec
220 call print_message( msgbuf, standardmessageunit,
221 & SQUEEZE_RIGHT , mythid)
222 write(msgbuf,'(a,a,d22.15)')
223 & ' global cost function value',
224 & ' ( CTD temp. ) = ',fcthread_ctdt
225 call print_message( msgbuf, standardmessageunit,
226 & SQUEEZE_RIGHT , mythid)
227 write(msgbuf,'(a)') ' '
228 call print_message( msgbuf, standardmessageunit,
229 & SQUEEZE_RIGHT , mythid)
230 #endif
231
232 enddo
233 c-- End of second loop over records.
234
235 #endif
236
237 return
238 end

  ViewVC Help
Powered by ViewVC 1.1.22