/[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.1 - (hide annotations) (download)
Wed May 22 19:39:51 2013 UTC (13 years, 3 months ago) by jmc
Branch: MAIN
unfinished set-up to test pkg/cheapaml with thsice

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

  ViewVC Help
Powered by ViewVC 1.1.22