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

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

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


Revision 1.9 - (show annotations) (download)
Tue Oct 9 00:02:51 2007 UTC (16 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59k, checkpoint59j, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.8: +7 -6 lines
add missing cvs $Header:$ or $Name:$

1 C $Header: $
2 C $Name: $
3
4 #include "COST_CPPOPTIONS.h"
5
6
7 subroutine cost_theta( myiter, mytime, mythid )
8
9 c ==================================================================
10 c SUBROUTINE cost_theta
11 c ==================================================================
12 c
13 c o Evaluate cost function contribution of temperature.
14 c
15 c started: Christian Eckert eckert@mit.edu 30-Jun-1999
16 c
17 c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
18 c
19 c - Restructured the code in order to create a package
20 c for the MITgcmUV.
21 c
22 c changed: Patrick Heimbach heimbach@mit.edu 27-May-2000
23 c
24 c - set ladinit to .true. to initialise adtbar file
25 c
26 c ==================================================================
27 c SUBROUTINE cost_theta
28 c ==================================================================
29
30 implicit none
31
32 c == global variables ==
33
34 #include "EEPARAMS.h"
35 #include "SIZE.h"
36 #include "PARAMS.h"
37 #include "GRID.h"
38 #include "DYNVARS.h"
39
40 #include "cal.h"
41 #include "ecco_cost.h"
42 #include "ctrl.h"
43 #include "ctrl_dummy.h"
44 #include "optim.h"
45
46 c == routine arguments ==
47
48 integer myiter
49 _RL mytime
50 integer mythid
51
52 c == local variables ==
53
54 integer bi,bj
55 integer i,j,k
56 integer itlo,ithi
57 integer jtlo,jthi
58 integer jmin,jmax
59 integer imin,imax
60 integer irec, irectmp
61 integer levmon
62 integer levoff
63 integer iltheta
64
65 _RL fctile
66 _RL fcthread
67
68 _RL cmask (1-olx:snx+olx,1-oly:sny+oly)
69 _RL spval
70 _RL spmax
71
72 character*(80) fnametheta
73
74 logical doglobalread
75 logical ladinit
76
77 character*(MAX_LEN_MBUF) msgbuf
78 #ifdef GENERIC_BAR_MONTH
79 integer mrec, nyears, iyear
80 #endif
81
82 _RL diagnosfld3d(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
83
84 c == external functions ==
85
86 integer ilnblnk
87 external ilnblnk
88
89 c == end of interface ==
90
91 jtlo = mybylo(mythid)
92 jthi = mybyhi(mythid)
93 itlo = mybxlo(mythid)
94 ithi = mybxhi(mythid)
95 jmin = 1
96 jmax = sny
97 imin = 1
98 imax = snx
99
100 spval = -1.8
101 spmax = 40.
102
103 c-- Read tiled data.
104 doglobalread = .false.
105 ladinit = .false.
106
107 #ifdef ALLOW_THETA_COST_CONTRIBUTION
108
109 if (optimcycle .ge. 0) then
110 iltheta = ilnblnk( tbarfile )
111 write(fnametheta(1:80),'(2a,i10.10)')
112 & tbarfile(1:iltheta),'.',optimcycle
113 endif
114
115 fcthread = 0. _d 0
116
117 #ifdef GENERIC_BAR_MONTH
118 c-- Loop over month
119 do irec = 1,min(nmonsrec,12)
120 nyears=int((nmonsrec-irec)/12)+1
121 do iyear=1,nyears
122 mrec=irec+(iyear-1)*12
123 irectmp=mrec
124 c-- Read time averages and the monthly mean data.
125 call active_read_xyz( fnametheta, tbar, mrec,
126 & doglobalread, ladinit,
127 & optimcycle, mythid,
128 & xx_tbar_mean_dummy )
129 do bj = jtlo,jthi
130 do bi = itlo,ithi
131 do k = 1,nr
132 do j = jmin,jmax
133 do i = imin,imax
134 if(iyear.eq.1) then
135 tbar_gen(i,j,k,bi,bj) =tbar(i,j,k,bi,bj)
136 elseif(iyear.eq.nyears) then
137 tbar(i,j,k,bi,bj) =(tbar_gen(i,j,k,bi,bj)
138 $ +tbar(i,j,k,bi,bj))/float(nyears)
139 else
140 tbar_gen(i,j,k,bi,bj) =tbar_gen(i,j,k,bi,bj)
141 $ +tbar(i,j,k,bi,bj)
142 endif
143 enddo
144 enddo
145 enddo
146 enddo
147 enddo
148 enddo
149 #else
150 c-- Loop over records.
151 do irec = 1,nmonsrec
152
153 irectmp = irec
154 c-- Read time averages and the monthly mean data.
155 call active_read_xyz( fnametheta, tbar, irec,
156 & doglobalread, ladinit,
157 & optimcycle, mythid,
158 & xx_tbar_mean_dummy )
159 #endif
160 c-- Determine the month to be read.
161 levoff = mod(modelstartdate(1)/100,100)
162 levmon = (irectmp-1) + levoff
163 levmon = mod(levmon-1,12)+1
164
165 call mdsreadfield( tdatfile, cost_iprec, cost_yftype,
166 & nr, tdat, levmon, mythid)
167
168 do bj = jtlo,jthi
169 do bi = itlo,ithi
170
171 c-- Loop over the model layers
172 fctile = 0. _d 0
173 do k = 1,nr
174
175 c-- Determine the mask on weights
176 do j = jmin,jmax
177 do i = imin,imax
178 cmask(i,j) = cosphi(i,j,bi,bj)
179 if (tdat(i,j,k,bi,bj) .eq. 0.) then
180 cmask(i,j) = 0. _d 0
181 else if (tdat(i,j,k,bi,bj) .lt. spval) then
182 cmask(i,j) = 0. _d 0
183 else if (tdat(i,j,k,bi,bj) .gt. spmax) then
184 cmask(i,j) = 0. _d 0
185 endif
186 enddo
187 enddo
188
189 c-- Compute model data misfit and cost function term for
190 c the temperature field.
191 do j = jmin,jmax
192 do i = imin,imax
193 if ( _hFacC(i,j,k,bi,bj) .ne. 0. ) then
194 fctile = fctile +
195 & (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
196 & (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
197 & (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
198 if ( wthetaLev(i,j,k,bi,bj)*cmask(i,j) .ne. 0. )
199 & num_temp(bi,bj) = num_temp(bi,bj) + 1. _d 0
200 diagnosfld3d(i,j,k,bi,bj) =
201 & (wthetaLev(i,j,k,bi,bj)*cmask(i,j)*
202 & (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj))*
203 & (tbar(i,j,k,bi,bj) - tdat(i,j,k,bi,bj)) )
204 else
205 diagnosfld3d(i,j,k,bi,bj) = 0.
206 endif
207 enddo
208 enddo
209
210 enddo
211 c-- End of loop over layers.
212
213 call mdswritefield( 'DiagnosCost_ClimTheta',
214 & writeBinaryPrec, globalfiles, 'RL', Nr,
215 & diagnosfld3d, irec, optimcycle, mythid )
216
217 fcthread = fcthread + fctile
218 objf_temp(bi,bj) = objf_temp(bi,bj) + fctile
219
220 enddo
221 enddo
222
223 enddo
224 c-- End of loop over records.
225
226 #endif
227
228 return
229 end
230

  ViewVC Help
Powered by ViewVC 1.1.22