/[MITgcm]/MITgcm/pkg/seaice/cost_ice_test.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/cost_ice_test.F

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


Revision 1.14 - (hide annotations) (download)
Wed Dec 3 16:47:45 2014 UTC (9 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, HEAD
Changes since 1.13: +4 -4 lines
replace a few print *, statements by write(standardMessageUnit,*),
because on some systems this means different output streams

1 mlosch 1.14 C $Header: /u/gcmpack/MITgcm/pkg/seaice/cost_ice_test.F,v 1.13 2014/04/28 15:09:27 mlosch Exp $
2 jmc 1.4 C $Name: $
3 heimbach 1.1
4     #include "SEAICE_OPTIONS.h"
5    
6     subroutine cost_ice_test( mytime, myiter, mythid )
7    
8     c ==================================================================
9     c SUBROUTINE cost_ice_test
10     c ==================================================================
11     c
12 dimitri 1.6 c o Compute sea-ice cost function. The following options can be
13     c selected with data.seaice (SEAICE_PARM02) variable cost_ice_flag
14 heimbach 1.1 c
15     c cost_ice_flag = 1
16     c - compute mean sea-ice volume
17     c costIceStart < mytime < costIceEnd
18     c
19     c cost_ice_flag = 2
20     c - compute mean sea-ice area
21     c costIceStart < mytime < costIceEnd
22     c
23     c cost_ice_flag = 3
24     c - heat content of top level plus latent heat of sea-ice
25     c costIceStart < mytime < costIceEnd
26     c
27     c cost_ice_flag = 4
28     c - heat content of top level
29     c costIceStart < mytime < costIceEnd
30     c
31     c cost_ice_flag = 5
32     c - heat content of top level plus sea-ice plus latent heat of snow
33     c costIceStart < mytime < costIceEnd
34     c
35     c cost_ice_flag = 6
36     c - quadratic cost function measuring difference between pkg/seaice
37     c AREA variable and simulated sea-ice measurements at every time
38     c step.
39     c
40     c ==================================================================
41     c
42     c started: menemenlis@jpl.nasa.gov 26-Feb-2003
43     c
44     c ==================================================================
45     c SUBROUTINE cost_ice_test
46     c ==================================================================
47    
48     implicit none
49    
50     c == global variables ==
51     #ifdef ALLOW_COST_ICE
52     #include "EEPARAMS.h"
53     #include "SIZE.h"
54     #include "GRID.h"
55     #include "PARAMS.h"
56 jmc 1.11 #include "SEAICE_SIZE.h"
57 heimbach 1.1 #include "SEAICE_COST.h"
58     #include "SEAICE.h"
59     #include "DYNVARS.h"
60     #include "cost.h"
61     #endif /* ALLOW_COST_ICE */
62    
63     c == routine arguments ==
64    
65     _RL mytime
66     integer myiter
67     integer mythid
68    
69 dimitri 1.5 #ifdef ALLOW_COST_ICE
70 heimbach 1.1
71     c == local variables ==
72    
73     c msgBuf - Informational/error message buffer
74     CHARACTER*(MAX_LEN_MBUF) msgBuf
75 mlosch 1.9 integer bi,bj,i,j,kSrf
76 heimbach 1.1 _RL tempVar
77    
78     c == external functions ==
79    
80     integer ilnblnk
81     external ilnblnk
82    
83     c == end of interface ==
84    
85     if ( myTime .GT. (endTime - lastinterval) ) then
86 jmc 1.10 tempVar = 1. _d 0/
87     & ( ( 1. _d 0 + min(endTime-startTime,lastinterval) )
88 jmc 1.4 & / deltaTClock )
89 heimbach 1.1
90 mlosch 1.9 kSrf = 1
91 heimbach 1.1 cph(
92 mlosch 1.14 write(standardMessageUnit,*) 'ph-ice B ', myiter,
93     & theta(4,4,kSrf,1,1), area(4,4,1,1), heff(4,4,1,1)
94 heimbach 1.1 cph)
95     if ( cost_ice_flag .eq. 1 ) then
96     c sea-ice volume
97     do bj=myByLo(myThid),myByHi(myThid)
98     do bi=myBxLo(myThid),myBxHi(myThid)
99     do j = 1,sny
100     do i = 1,snx
101     objf_ice(bi,bj) = objf_ice(bi,bj) +
102 mlosch 1.8 & tempVar * rA(i,j,bi,bj) * HEFF(i,j,bi,bj)
103 heimbach 1.1 enddo
104     enddo
105     enddo
106     enddo
107    
108     elseif ( cost_ice_flag .eq. 2 ) then
109     c sea-ice area
110     do bj=myByLo(myThid),myByHi(myThid)
111     do bi=myBxLo(myThid),myBxHi(myThid)
112     do j = 1,sny
113     do i = 1,snx
114     objf_ice(bi,bj) = objf_ice(bi,bj) +
115 mlosch 1.8 & tempVar * rA(i,j,bi,bj) * AREA(i,j,bi,bj)
116 heimbach 1.1 enddo
117     enddo
118     enddo
119     enddo
120    
121     c heat content of top level:
122     c theta * delZ * (sea water heat capacity = 3996 J/kg/K)
123     c * (density of sea-water = 1026 kg/m^3)
124     c
125     c heat content of sea-ice:
126     c tice * heff * (sea ice heat capacity = 2090 J/kg/K)
127     c * (density of sea-ice = 910 kg/m^3)
128     c
129     c note: to remove mass contribution to heat content,
130     c which is not properly accounted for by volume converving
131     c ocean model, theta and tice are referenced to freezing
132     c temperature of sea-ice, here -1.96 deg C
133     c
134     c latent heat content of sea-ice:
135     c - heff * (latent heat of fusion = 334000 J/kg)
136     c * (density of sea-ice = 910 kg/m^3)
137     c
138     c latent heat content of snow:
139     c - hsnow * (latent heat of fusion = 334000 J/kg)
140     c * (density of snow = 330 kg/m^3)
141    
142     elseif ( cost_ice_flag .eq. 3 ) then
143     c heat content of top level plus latent heat of sea-ice
144     do bj=myByLo(myThid),myByHi(myThid)
145     do bi=myBxLo(myThid),myBxHi(myThid)
146     do j = 1,sny
147     do i = 1,snx
148     objf_ice(bi,bj) = objf_ice(bi,bj) +
149     & tempVar * rA(i,j,bi,bj) * (
150 jmc 1.10 & (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
151     & drF(1) * 3996. _d 0 * 1026. _d 0 -
152     & HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 )
153 heimbach 1.1 enddo
154     enddo
155     enddo
156     enddo
157    
158     elseif ( cost_ice_flag .eq. 4 ) then
159     c heat content of top level
160     do bj=myByLo(myThid),myByHi(myThid)
161     do bi=myBxLo(myThid),myBxHi(myThid)
162     do j = 1,sny
163     do i = 1,snx
164     objf_ice(bi,bj) = objf_ice(bi,bj) +
165     & tempVar * rA(i,j,bi,bj) * (
166 jmc 1.10 & (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
167     & drF(1) * 3996. _d 0 * 1026. _d 0 )
168 heimbach 1.1 enddo
169     enddo
170     enddo
171     enddo
172    
173     elseif ( cost_ice_flag .eq. 5 ) then
174     c heat content of top level plus sea-ice plus latent heat of snow
175     do bj=myByLo(myThid),myByHi(myThid)
176     do bi=myBxLo(myThid),myBxHi(myThid)
177     do j = 1,sny
178     do i = 1,snx
179     objf_ice(bi,bj) = objf_ice(bi,bj) +
180     & tempVar * rA(i,j,bi,bj) * (
181 jmc 1.10 & (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
182     & drF(1) * 3996. _d 0 * 1026. _d 0 +
183 mlosch 1.13 & (TICES(i,j,1,bi,bj) - 273.15 _d 0 + 1.96 _d 0 ) *
184 jmc 1.10 & HEFF(i,j,bi,bj) * 2090. _d 0 * 910. _d 0 -
185     & HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 -
186     & HSNOW(i,j,bi,bj) * 334000. _d 0 * 330. _d 0 )
187 heimbach 1.1 enddo
188     enddo
189     enddo
190     enddo
191    
192     elseif ( cost_ice_flag .eq. 6 ) then
193     c Qadratic cost function measuring difference between pkg/seaice
194     c AREA variable and simulated sea-ice measurements at every time
195     c step. For time being no measurements are read-in. It is
196     c assumed that measurements are AREA=0.5 at all times everywhere.
197     do bj=myByLo(myThid),myByHi(myThid)
198     do bi=myBxLo(myThid),myBxHi(myThid)
199     do j = 1,sny
200     do i = 1,snx
201     objf_ice(bi,bj) = objf_ice(bi,bj) +
202 jmc 1.10 & ( AREA(i,j,bi,bj) - 0.5 _d 0 ) *
203     & ( AREA(i,j,bi,bj) - 0.5 _d 0 )
204 heimbach 1.1 enddo
205     enddo
206     enddo
207     enddo
208    
209 heimbach 1.12 elseif ( cost_ice_flag .eq. 7 ) then
210    
211     do bj=myByLo(myThid),myByHi(myThid)
212     do bi=myBxLo(myThid),myBxHi(myThid)
213     do j = 1,sny
214     do i = 1,snx
215     objf_ice(bi,bj) = objf_ice(bi,bj) +
216     & UICE(i,j,bi,bj) * UICE(i,j,bi,bj) +
217     & VICE(i,j,bi,bj) * VICE(i,j,bi,bj)
218    
219     enddo
220     enddo
221     enddo
222     enddo
223    
224 heimbach 1.1 else
225     WRITE(msgBuf,'(A)')
226     & 'COST_ICE: invalid cost_ice_flag'
227     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
228     & SQUEEZE_RIGHT , myThid )
229     STOP 'ABNORMAL END: S/R COST_ICE'
230     endif
231     endif
232    
233 heimbach 1.2 cph(
234 mlosch 1.14 write(standardMessageUnit,*) 'ph-ice C ', myiter, objf_ice(1,1)
235 heimbach 1.2 cph)
236    
237 heimbach 1.1 #endif /* ALLOW_COST_ICE */
238    
239 jmc 1.10 return
240 heimbach 1.1 end

  ViewVC Help
Powered by ViewVC 1.1.22