/[MITgcm]/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml_seaice.F
ViewVC logotype

Annotation of /MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml_seaice.F

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

Revision 1.2 - (hide annotations) (download)
Sat May 25 18:19:06 2013 UTC (13 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64h
Changes since 1.1: +56 -21 lines 
- consistent treatment of uWind,vWind location: assume everywhere in
  pkg/cheapaml that they are on C-grid (@ uVel,vVel location)
  (Note: was already the case for Tair,Qair advection);
- fill in array "Qsw" (for short-wave heating);
- add diagnostic for precip ;
1 jmc 1.2 C $Header: /u/gcmpack/MITgcm_contrib/verification_other/offline_cheapaml/code/cheapaml_seaice.F,v 1.1 2013/05/22 19:39:51 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CHEAPAML_OPTIONS.h"
5     #ifdef ALLOW_THSICE
6     # include "THSICE_OPTIONS.h"
7     #endif
8 jmc 1.2 #ifdef ALLOW_SEAICE
9     # include "SEAICE_OPTIONS.h"
10     #endif
11 jmc 1.1
12     CBOP
13     C !ROUTINE: CHEAPAML_SEAICE
14     C !INTERFACE:
15     SUBROUTINE CHEAPAML_SEAICE(
16 jmc 1.2 I swDown, lwDown, uWind, vWind, LVapor,
17 jmc 1.1 O fsha, flha, evp, xolw, ssqt, q100, cdq,
18     O Tsurf, iceFrac, sw2oce,
19     I bi, bj, myTime, myIter, myThid )
20     C !DESCRIPTION: \bv
21     C *==========================================================*
22     C | S/R CHEAPAML_SEAICE
23     C | o Compute fluxes over seaice by calling seaice routine
24     C | to solve for surface temperature.
25     C *==========================================================*
26     C \ev
27    
28     C !USES:
29     IMPLICIT NONE
30     C == Global variables ===
31     #include "SIZE.h"
32     #include "EEPARAMS.h"
33     #include "PARAMS.h"
34     #ifdef ALLOW_THSICE
35     #include "THSICE_PARAMS.h"
36     #include "THSICE_SIZE.h"
37     #include "THSICE_VARS.h"
38 jmc 1.2 #endif
39     #ifdef ALLOW_SEAICE
40     # include "SEAICE_SIZE.h"
41     # include "SEAICE.h"
42 jmc 1.1 #endif
43    
44     INTEGER siLo, siHi, sjLo, sjHi
45     PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
46     PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
47    
48     C !INPUT PARAMETERS:
49     C == Routine Arguments ==
50     C swDown :: incoming short-wave radiation (+=dw) [W/m2]
51     C lwDown :: incoming long-wave radiation (+=dw) [W/m2]
52     C uRelWind :: relative wind speed, u-component [m/s]
53     C vRelWind :: relative wind speed, v-component [m/s]
54     C LVapor :: latent heat of vaporisation
55     C bi, bj :: tile indices
56     C myIter :: current iteration number
57     C myTime :: current time in simulation
58     C myThid :: my Thread Id number
59     _RL swDown(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60     _RL lwDown(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 jmc 1.2 _RL uWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62     _RL vWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63 jmc 1.1 _RL LVapor
64     _RL myTime
65     INTEGER bi, bj, myIter, myThid
66    
67     C !OUTPUT PARAMETERS:
68     C fsha :: sensible heat-flux over seaice (+=up) [W/m2]
69     C flha :: latent heat-flux over seaice (+=up) [W/m2]
70     C evp :: evaporation over seaice (+=up) [kg/m2/s]
71     C xolw :: upward long-wave over seaice (+=up) [W/m2]
72     C ssqt ::
73     C q100 ::
74     C cdq ::
75     C Tsurf :: updated seaice/snow surface temperature [deg.C]
76     C iceFrac :: ice fraction [0-1]
77     C sw2oce :: short-wave over seaice into the ocean (+=dw) [W/m2]
78     _RL fsha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79     _RL flha(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80     _RL evp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81     _RL xolw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82     _RL ssqt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83     _RL q100(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84     _RL cdq (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85     _RL Tsurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86     _RL iceFrac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87     _RS sw2oce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88     c _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
89    
90     #ifdef ALLOW_THSICE
91     C !LOCAL VARIABLES:
92     C == Local variables ==
93 jmc 1.2 C uRelWind :: relative wind speed, u-component [m/s], (C-grid)
94     C vRelWind :: relative wind speed, v-component [m/s], (C-grid)
95     C windSq :: relative wind speed squared (grid-cell center)
96 jmc 1.1 INTEGER i, j
97     INTEGER iceOrNot
98     INTEGER iMin, iMax
99     INTEGER jMin, jMax
100     _RL LatentHeat
101     _RL icFrac, opFrac
102     _RL netSW (1:sNx,1:sNy)
103     _RL sFlx (1:sNx,1:sNy,0:2)
104     c _RL tFrzOce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL dTsurf(1:sNx,1:sNy)
106    
107 jmc 1.2 _RL uRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL vRelWind(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     _RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110     _RL cdu, dumArg(4)
111 jmc 1.1 _RL fsha0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112     _RL evp_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113     _RL xolw0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114     c _RL ssqt0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115     c _RL q10_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116     c _RL cdq_0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117     _RL dShdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118     _RL dEvdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119     _RL dLwdTs(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120     CEOP
121    
122     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
123    
124     iMin = 1
125     iMax = sNx
126     jMin = 1
127     jMax = sNy
128     LatentHeat = Lfresh + LVapor
129    
130     c DO bj=myByLo(myThid),myByHi(myThid)
131     c DO bi=myBxLo(myThid),myBxHi(myThid)
132    
133     CALL THSICE_GET_OCEAN(
134     I bi, bj, myTime, myIter, myThid )
135    
136 jmc 1.2 DO j = 1-OLy, sNy+OLy
137     DO i = 1-OLx, sNx+OLx
138     uRelWind(i,j) = uWind(i,j)
139     vRelWind(i,j) = vWind(i,j)
140     ENDDO
141     ENDDO
142     #ifdef ALLOW_SEAICE
143     IF ( useSEAICE ) THEN
144     DO j = 1-OLy, sNy+OLy
145     DO i = 1-OLx, sNx+OLx
146     uRelWind(i,j) = uRelWind(i,j)-uIce(i,j,bi,bj)
147     vRelWind(i,j) = vRelWind(i,j)-vIce(i,j,bi,bj)
148     ENDDO
149     ENDDO
150     ENDIF
151     #endif /* ALLOW_SEAICE */
152     DO j = jMin,jMax
153     DO i = iMin,iMax
154     windSq(i,j) = ( uRelWind( i ,j)*uRelWind( i ,j)
155     & + uRelWind(i+1,j)*uRelWind(i+1,j)
156     & + vRelWind(i, j )*vRelWind(i, j )
157     & + vRelWind(i,j+1)*vRelWind(i,j+1)
158     & )*0.5 _d 0
159     ENDDO
160     ENDDO
161    
162 jmc 1.1 C 1) compute albedo ; compute netSW
163     CALL THSICE_ALBEDO(
164     I bi, bj, siLo, siHi, sjLo, sjHi,
165     I iMin,iMax, jMin,jMax,
166     I iceMask(siLo,sjLo,bi,bj), iceHeight(siLo,sjLo,bi,bj),
167     I snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
168     I snowAge(siLo,sjLo,bi,bj),
169     O siceAlb(siLo,sjLo,bi,bj), icAlbNIR(siLo,sjLo,bi,bj),
170     I myTime, myIter, myThid )
171    
172     DO j = jMin, jMax
173     DO i = iMin, iMax
174     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
175     C- surface net SW flux:
176     netSW(i,j) = swDown(i,j)
177     & *(1. _d 0 - siceAlb(i,j,bi,bj))
178     ELSE
179     netSW(i,j) = swDown(i,j)
180     ENDIF
181     ENDDO
182     ENDDO
183    
184    
185     C 2) compute other flx over seaice, over melting surf
186     C 3) compute other flx over seaice & derivative vs Tsurf, using previous Tsurf
187     DO j = jMin, jMax
188     DO i = iMin, iMax
189    
190     IF ( snowHeight(i,j,bi,bj).GT.3. _d -1 ) THEN
191     iceornot=2
192     ELSE
193     iceornot=1
194     ENDIF
195     Tsurf(i,j) = 0.
196     CALL CHEAPAML_COARE3_FLUX(
197 jmc 1.2 I i, j, bi, bj, iceOrNot,
198     I Tsurf, windSq,
199     O fsha0(i,j), flha(i,j), evp_0(i,j), xolw0(i,j),
200     O ssqt(i,j), q100(i,j), cdq(i,j), cdu,
201     O dumArg(1), dumArg(2), dumArg(3), dumArg(4),
202     I myIter, myThid )
203 jmc 1.1 sFlx(i,j,0) = lwDown(i,j)- xolw0(i,j)
204     & - fsha0(i,j) - evp_0(i,j)*LatentHeat
205    
206     Tsurf(i,j) = Tsrf(i,j,bi,bj)
207     CALL CHEAPAML_COARE3_FLUX(
208 jmc 1.2 I i, j, bi, bj, iceOrNot,
209     I Tsurf, windSq,
210     O fsha(i,j), flha(i,j), evp(i,j), xolw(i,j),
211     O ssqt(i,j), q100(i,j), cdq(i,j), cdu,
212     O dShdTs(i,j), dEvdTs(i,j), dLwdTs(i,j), dumArg(4),
213     I myIter, myThid )
214 jmc 1.1 sFlx(i,j,1) = lwDown(i,j)- xolw(i,j)
215     & - fsha(i,j) - evp(i,j)*LatentHeat
216     sFlx(i,j,2) = -dLwdTs(i,j)
217     & - dShdTs(i,j) - dEvdTs(i,j)*LatentHeat
218     ENDDO
219     ENDDO
220    
221     C 4) solve for surf & seaice temp
222     C-- needs to fill in snowPrc, ( & prcAtm ? )
223     C-- note: this S/R assumes No overlap
224     CALL THSICE_IMPL_TEMP(
225     I netSW, sFlx,
226     O dTsurf,
227     I bi, bj, myTime, myIter, myThid )
228    
229     C 5) update surf fluxes
230     DO j = jMin, jMax
231     DO i = iMin, iMax
232     iceFrac(i,j) = iceMask(i,j,bi,bj)
233     sw2oce (i,j) = icFlxSW(i,j,bi,bj)
234     IF ( dTsurf(i,j) .GT. 999. ) THEN
235     c dTsurf(J)= tFreeze - Tsurf(J)
236     Tsurf(i,j)= 0.
237     fsha(i,j) = fsha0(i,j)
238     flha(i,j) = evp_0(i,j)*LatentHeat
239     evp(i,j) = evp_0(i,j)
240     xolw(i,j) = xolw0(i,j)
241     ELSE
242     Tsurf(i,j)= Tsurf(i,j)+ dTsurf(i,j)
243     fsha(i,j) = fsha(i,j) + dTsurf(i,j)*dShdTs(i,j)
244     evp(i,j) = evp(i,j) + dTsurf(i,j)*dEvdTs(i,j)
245     flha(i,j) = evp(i,j)*LatentHeat
246     xolw(i,j) = xolw(i,j) + dTsurf(i,j)*dLwdTs(i,j)
247     ENDIF
248     ENDDO
249     ENDDO
250    
251     DO j = jMin, jMax
252     DO i = iMin, iMax
253     c IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
254     icFrac = iceMask(i,j,bi,bj)
255     opFrac = 1. _d 0 - icFrac
256     C-- Update Fluxes :
257     icFlxAtm(i,j,bi,bj) = netSW(i,j)
258     & + lwDown(i,j)- xolw(i,j)
259     & - fsha(i,j) - evp(i,j)*LVapor
260     icFrwAtm(i,j,bi,bj) = evp(i,j)
261     c ENDIF
262     ENDDO
263     ENDDO
264    
265     c ENDDO
266     c ENDDO
267    
268     #endif /* ALLOW_THSICE */
269     RETURN
270     END
  ViewVC Help
Powered by ViewVC 1.1.22