/[MITgcm]/MITgcm_contrib/SOSE/code_ad/seaice_cost_areasst.F
ViewVC logotype

Annotation of /MITgcm_contrib/SOSE/code_ad/seaice_cost_areasst.F

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


Revision 1.1 - (hide annotations) (download)
Fri Apr 23 19:55:12 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD
original files

1 mmazloff 1.1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_cost_areasst.F,v 1.2 2008/10/14 17:41:40 jmc Exp $
2     C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5     #ifdef ALLOW_COST
6     #include "COST_CPPOPTIONS.h"
7     #endif
8    
9    
10     subroutine seaice_cost_areasst(
11     & nnzbar, localbarfile, localbar, xx_localbar_mean_dummy,
12     & nnzobs, localobsfile, localobs, mult_local,
13     & nrecloc, localstartdate, localperiod,
14     & localmask, localweight,
15     & spminloc, spmaxloc, spzeroloc,
16     & objf_local, num_local,
17     & myiter, mytime, mythid )
18    
19     c ==================================================================
20     c SUBROUTINE seaice_cost_areasst
21     c ==================================================================
22     c
23     c o Based on cost_generic
24     c o in case where there is observed sea-ice we not only constrain
25     c model(area)=obs(area) but also model(sst)@freezing point
26     c
27     c ==================================================================
28     c SUBROUTINE seaice_cost_areasst
29     c ==================================================================
30    
31     implicit none
32    
33     c == global variables ==
34    
35     #include "EEPARAMS.h"
36     #include "SIZE.h"
37     #include "PARAMS.h"
38     #ifdef ALLOW_CAL
39     # include "cal.h"
40     #endif
41     #include "SEAICE.h"
42     #include "SEAICE_PARAMS.h"
43     #include "SEAICE_COST.h"
44     #ifdef ALLOW_COST
45     # ifdef ALLOW_ECCO
46     # include "ecco_cost.h"
47     # endif
48     # include "optim.h"
49     #endif
50     #ifdef ALLOW_CTRL
51     # include "ctrl.h"
52     # include "ctrl_dummy.h"
53     #endif
54    
55     c == routine arguments ==
56    
57     integer nnzbar
58     integer nnzobs
59     integer nrecloc
60     integer myiter
61     integer mythid
62     integer localstartdate(4)
63    
64     _RL localbarT (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
65     _RL localbar (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
66     _RL localobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
67     _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
68     _RL localmask (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
69     _RL xx_localbar_mean_dummy
70     _RL mult_local
71     _RL mytime
72     _RL localperiod
73     _RL spminloc
74     _RL spmaxloc
75     _RL spzeroloc
76     _RL objf_local(nsx,nsy)
77     _RL num_local(nsx,nsy)
78    
79     character*(MAX_LEN_FNAM) localbarfile
80     character*(MAX_LEN_FNAM) localobsfile
81    
82     #ifdef ALLOW_ECCO
83     #ifdef ALLOW_COST
84     c == local variables ==
85    
86     integer bi,bj
87     integer i,j,k
88     integer itlo,ithi
89     integer jtlo,jthi
90     integer jmin,jmax
91     integer imin,imax
92     integer irec
93     integer il
94     integer localrec
95     integer obsrec
96    
97     logical doglobalread
98     logical ladinit
99    
100     _RL spval
101     parameter (spval = -9999. )
102     _RL localwww
103     _RL localcost
104     _RL junk
105    
106     _RL cmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs)
107    
108     character*(128) fname1, fname2
109     character*(MAX_LEN_MBUF) msgbuf
110    
111     cnew(
112     _RL daytime
113     _RL diffsecs
114     integer dayiter
115     integer daydate(4)
116     integer difftime(4)
117     integer middate(4)
118     integer yday, ymod
119     integer md, dd, sd, ld, wd
120     integer mody, modm
121     integer beginmodel, beginlocal
122     logical exst
123     cnew)
124    
125     c == external functions ==
126    
127     integer ilnblnk
128     external ilnblnk
129    
130     c == end of interface ==
131    
132     jtlo = mybylo(mythid)
133     jthi = mybyhi(mythid)
134     itlo = mybxlo(mythid)
135     ithi = mybxhi(mythid)
136     jmin = 1
137     jmax = sny
138     imin = 1
139     imax = snx
140    
141     c-- Initialise local variables.
142    
143     localwww = 0. _d 0
144    
145     do bj = jtlo,jthi
146     do bi = itlo,ithi
147     objf_local(bi,bj) = 0. _d 0
148     num_local(bi,bj) = 0. _d 0
149     enddo
150     enddo
151    
152     c-- First, read tiled data.
153     doglobalread = .false.
154     ladinit = .false.
155    
156     write(fname1(1:128),'(80a)') ' '
157     il=ilnblnk( localbarfile )
158     write(fname1(1:128),'(2a,i10.10)')
159     & localbarfile(1:il),'.',optimcycle
160    
161     cph if ( .NOT. ( mult_local.EQ.0. .OR. localobsfile.EQ.' ' ) ) then
162     if ( .NOT. ( localobsfile.EQ.' ' ) ) then
163    
164     c-- Loop over records for the second time.
165     do irec = 1, nrecloc
166    
167     if ( nnzbar .EQ. 1 ) then
168     call active_read_xy( fname1, localbar, irec, doglobalread,
169     & ladinit, optimcycle, mythid,
170     & xx_localbar_mean_dummy )
171     else
172     call active_read_xyz( fname1, localbar, irec, doglobalread,
173     & ladinit, optimcycle, mythid,
174     & xx_localbar_mean_dummy )
175     endif
176    
177     cnew(
178     if ( localperiod .EQ. 86400. ) then
179     c-- assume daily fields
180     obsrec = irec
181     daytime = FLOAT(secondsperday*(irec-1))
182     dayiter = hoursperday*(irec-1)
183     call cal_getdate( dayiter, daytime, daydate, mythid )
184     call cal_convdate( daydate,yday,md,dd,sd,ld,wd,mythid )
185     ymod = localstartdate(1)/10000
186     if ( ymod .EQ. yday ) then
187     middate(1) = modelstartdate(1)
188     else
189     middate(1) = yday*10000+100+1
190     endif
191     middate(2) = 0
192     middate(3) = modelstartdate(3)
193     middate(4) = modelstartdate(4)
194     call cal_TimePassed( middate, daydate, difftime, mythid )
195     call cal_ToSeconds( difftime, diffsecs, mythid )
196     localrec = int(diffsecs/localperiod) + 1
197     else
198     c-- assume monthly fields
199     beginlocal = localstartdate(1)/10000
200     beginmodel = modelstartdate(1)/10000
201     obsrec =
202     & ( beginmodel - beginlocal )*nmonthyear
203     & + ( mod(modelstartdate(1)/100,100)
204     & -mod(localstartdate(1)/100,100) )
205     & + irec
206     mody = modelstartdate(1)/10000
207     modm = modelstartdate(1)/100 - mody*100
208     yday = mody + INT((modm-1+irec-1)/12)
209     localrec = 1 + MOD(modm-1+irec-1,12)
210     endif
211    
212     il=ilnblnk(localobsfile)
213     write(fname2(1:128),'(2a,i4)')
214     & localobsfile(1:il), '_', yday
215     inquire( file=fname2, exist=exst )
216     if (.NOT. exst) then
217     write(fname2(1:128),'(a)') localobsfile(1:il)
218     localrec = obsrec
219     endif
220     if ( localrec .GT. 0 ) then
221     call mdsreadfield( fname2, cost_iprec, cost_yftype, nnzobs,
222     & localobs, localrec, mythid )
223     else
224     do bj = jtlo,jthi
225     do bi = itlo,ithi
226     do k = 1,nnzobs
227     do j = jmin,jmax
228     do i = imin,imax
229     localobs(i,j,bi,bj) = spval
230     enddo
231     enddo
232     enddo
233     enddo
234     enddo
235     endif
236     cnew)
237    
238     do bj = jtlo,jthi
239     do bi = itlo,ithi
240    
241     localcost = 0. _d 0
242    
243     c-- Determine the mask on weights
244     do k = 1,nnzobs
245     do j = jmin,jmax
246     do i = imin,imax
247     cmask(i,j,k) = cosphi(i,j,bi,bj)*localmask(i,j,k,bi,bj)
248     if ( localobs(i,j,bi,bj) .lt. spminloc .or.
249     & localobs(i,j,bi,bj) .gt. spmaxloc .or.
250     & localobs(i,j,bi,bj) .eq. spzeroloc ) then
251     cmask(i,j,k) = 0. _d 0
252     endif
253     enddo
254     enddo
255     enddo
256     c--
257     do k = 1,nnzobs
258     do j = jmin,jmax
259     do i = imin,imax
260     localwww = localweight(i,j,bi,bj)*cmask(i,j,k)
261     CMM junk = ( localbar(i,j,k,bi,bj) -
262     CMM & localobs(i,j,k,bi,bj) )
263     junk = ( localbar(i,j,bi,bj) -
264     & SEAICE_freeze + 0.0001)
265     & *localobs(i,j,bi,bj)
266     localcost = localcost + junk*junk*localwww
267     if ( localwww .ne. 0. )
268     & num_local(bi,bj) = num_local(bi,bj) + 1. _d 0
269     enddo
270     enddo
271     enddo
272     objf_local(bi,bj) = objf_local(bi,bj) + localcost
273    
274     enddo
275     enddo
276    
277     enddo
278     c-- End of second loop over records.
279    
280     c-- End of mult_local or localobsfile
281     endif
282    
283     #endif /* ifdef ALLOW_COST */
284     #endif /* ifdef ALLOW_ECCO */
285    
286     end

  ViewVC Help
Powered by ViewVC 1.1.22