/[MITgcm]/MITgcm/verification/aim.5l_Equatorial_Channel/code/aim_surf_bc.F
ViewVC logotype

Contents of /MITgcm/verification/aim.5l_Equatorial_Channel/code/aim_surf_bc.F

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


Revision 1.8 - (show annotations) (download)
Thu Jan 21 00:13:14 2010 UTC (14 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint63, checkpoint62c, checkpoint62b, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.7: +77 -3 lines
update local version (more similar to std version)

1 C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/aim_surf_bc.F,v 1.16 2009/06/30 16:03:44 dfer Exp $
2 C $Name: $
3
4 #include "AIM_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: AIM_SURF_BC
8 C !INTERFACE:
9 SUBROUTINE AIM_SURF_BC(
10 U tYear,
11 O aim_sWght0, aim_sWght1,
12 I bi, bj, myTime, myIter, myThid )
13
14 C !DESCRIPTION: \bv
15 C *================================================================*
16 C | S/R AIM_SURF_BC
17 C | Set surface Boundary Conditions
18 C | for the atmospheric physics package
19 C *================================================================*
20 c | was part of S/R FORDATE in Franco Molteni SPEEDY code (ver23).
21 C | For now, surface BC are loaded from files (S/R AIM_FIELDS_LOAD)
22 C | and imposed (= surf. forcing).
23 C | In the future, will add
24 C | a land model and a coupling interface with an ocean GCM
25 C *================================================================*
26 C \ev
27
28 C !USES:
29 IMPLICIT NONE
30
31 C -------------- Global variables --------------
32 C-- size for MITgcm & Physics package :
33 #include "AIM_SIZE.h"
34
35 C-- MITgcm
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 c #include "DYNVARS.h"
39 #include "GRID.h"
40 c #include "SURFACE.h"
41
42 C-- Physics package
43 #include "AIM_PARAMS.h"
44 #include "AIM_FFIELDS.h"
45 c #include "AIM_GRID.h"
46 #include "com_forcon.h"
47 #include "com_forcing.h"
48 c #include "com_physvar.h"
49 #include "AIM_CO2.h"
50
51 C-- Coupled to the Ocean :
52 #ifdef COMPONENT_MODULE
53 #include "CPL_PARAMS.h"
54 #include "ATMCPL.h"
55 #endif
56
57 C !INPUT/OUTPUT PARAMETERS:
58 C == Routine arguments ==
59 C tYear :: Fraction into year
60 C aim_sWght0 :: weight for time interpolation of surface BC
61 C aim_sWght1 :: 0/1 = time period before/after the current time
62 C bi,bj :: Tile indices
63 C myTime :: Current time of simulation ( s )
64 C myIter :: Current iteration number in simulation
65 C myThid :: my Thread number Id.
66 _RL tYear
67 _RL aim_sWght0, aim_sWght1
68 INTEGER bi, bj
69 _RL myTime
70 INTEGER myIter, myThid
71 CEOP
72
73 #ifdef ALLOW_AIM
74 C !FUNCTIONS:
75 C !LOCAL VARIABLES:
76 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
77 C-- Local Variables originally (Speedy) in common bloc (com_forcing.h):
78 C-- COMMON /FORDAY/ Daily forcing fields (updated in FORDATE)
79 C oice1 :: sea ice fraction
80 C snow1 :: snow depth (mm water)
81 _RL oice1(NGP)
82 _RL snow1(NGP)
83 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
84 C == Local variables ==
85 C i,j,k,I2,k :: Loop counters
86 INTEGER i,j,I2,k, nm0
87 _RL t0prd, tNcyc, tmprd, dTprd
88 _RL SDEP1, IDEP2, SDEP2, SWWIL2, RSW, soilw_0, soilw_1
89 _RL RSD, alb_land, oceTfreez, ALBSEA1, ALPHA, CZEN, CZEN2
90 _RL RZEN, ZS, ZC, SJ, CJ, TMPA, TMPB, TMPL, hlim
91 c _RL DALB, alb_sea
92 #ifdef ALLOW_AIM_CO2
93 #ifdef ALLOW_DIAGNOSTICS
94 _RL pCO2scl
95 #endif
96 #endif /* ALLOW_AIM_CO2 */
97
98 C_EqCh: start
99 CHARACTER*(MAX_LEN_MBUF) suff
100 _RL xBump, yBump, dxBump, dyBump
101 xBump = xgOrigin + delX(1)*64.
102 yBump = ygOrigin + delY(1)*11.5
103 dxBump= delX(1)*12.
104 dyBump= delY(1)*6.
105 C_EqCh: Fix solar insolation with Sun directly overhead on Equator
106 tYear = 0.25 _d 0 - 10. _d 0/365. _d 0
107 C_EqCh: end
108
109 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
110 C- Set Land-sea mask (in [0,1]) from aim_landFr to fMask1:
111 DO j=1,sNy
112 DO i=1,sNx
113 I2 = i+(j-1)*sNx
114 fMask1(I2,1,myThid) = aim_landFr(i,j,bi,bj)
115 ENDDO
116 ENDDO
117
118 IF (aim_useFMsurfBC) THEN
119 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
120
121 C aim_surfForc_TimePeriod :: Length of forcing time period (e.g. 1 month)
122 C aim_surfForc_NppCycle :: Number of time period per Cycle (e.g. 12)
123 C aim_surfForc_TransRatio ::
124 C- define how fast the (linear) transition is from one month to the next
125 C = 1 -> linear between 2 midle month
126 C > TimePeriod/deltaT -> jump from one month to the next one
127
128 C-- Calculate weight for linear interpolation between 2 month centers
129 t0prd = myTime / aim_surfForc_TimePeriod
130 tNcyc = aim_surfForc_NppCycle
131 tmprd = t0prd - 0.5 _d 0 + tNcyc
132 tmprd = MOD(tmprd,tNcyc)
133 C- indices of previous month (nm0) and next month (nm1):
134 nm0 = 1 + INT(tmprd)
135 c nm1 = 1 + MOD(nm0,aim_surfForc_NppCycle)
136 C- interpolation weight:
137 dTprd = tmprd - (nm0 - 1)
138 aim_sWght1 = 0.5 _d 0+(dTprd-0.5 _d 0)*aim_surfForc_TransRatio
139 aim_sWght1 = MAX( 0. _d 0, MIN(1. _d 0, aim_sWght1) )
140 aim_sWght0 = 1. _d 0 - aim_sWght1
141
142 C-- Compute surface forcing at present time (linear Interp in time)
143 C using F.Molteni surface BC form ; fields needed are:
144 C 1. Sea Surface temperatures (in situ Temp. [K])
145 C 2. Land Surface temperatures (in situ Temp. [K])
146 C 3. Soil moisture (between 0-1)
147 C 4. Snow depth, Sea Ice : used to compute albedo (=> local arrays)
148 C 5. Albedo (between 0-1)
149
150 C- Surface Temperature:
151 DO j=1,sNy
152 DO i=1,sNx
153 I2 = i+(j-1)*sNx
154 sst1(I2,myThid) = aim_sWght0*aim_sst0(i,j,bi,bj)
155 & + aim_sWght1*aim_sst1(i,j,bi,bj)
156 stl1(I2,myThid) = aim_sWght0*aim_lst0(i,j,bi,bj)
157 & + aim_sWght1*aim_lst1(i,j,bi,bj)
158 ENDDO
159 ENDDO
160
161 C- Soil Water availability : (from F.M. INFORC S/R)
162 SDEP1 = 70. _d 0
163 IDEP2 = 3. _d 0
164 SDEP2 = IDEP2*SDEP1
165
166 SWWIL2= SDEP2*SWWIL
167 RSW = 1. _d 0/(SDEP1*SWCAP+SDEP2*(SWCAP-SWWIL))
168
169 DO j=1,sNy
170 DO i=1,sNx
171 I2 = i+(j-1)*sNx
172 soilw_0 = ( aim_sw10(i,j,bi,bj)
173 & +aim_veget(i,j,bi,bj)*
174 & MAX(IDEP2*aim_sw20(i,j,bi,bj)-SWWIL2, 0. _d 0)
175 & )*RSW
176 soilw_1 = ( aim_sw11(i,j,bi,bj)
177 & +aim_veget(i,j,bi,bj)*
178 & MAX(IDEP2*aim_sw21(i,j,bi,bj)-SWWIL2, 0. _d 0)
179 & )*RSW
180 soilw1(I2,myThid) = aim_sWght0*soilw_0
181 & + aim_sWght1*soilw_1
182 soilw1(I2,myThid) = MIN(1. _d 0, soilw1(I2,myThid) )
183 ENDDO
184 ENDDO
185
186 C- Set snow depth & sea-ice fraction :
187 DO j=1,sNy
188 DO i=1,sNx
189 I2 = i+(j-1)*sNx
190 snow1(I2) = aim_sWght0*aim_snw0(i,j,bi,bj)
191 & + aim_sWght1*aim_snw1(i,j,bi,bj)
192 oice1(I2) = aim_sWght0*aim_oic0(i,j,bi,bj)
193 & + aim_sWght1*aim_oic1(i,j,bi,bj)
194 ENDDO
195 ENDDO
196
197 IF (aim_splitSIOsFx) THEN
198 C- Split Ocean and Sea-Ice surf. temp. ; remove ice-fraction < 1 %
199 c oceTfreez = tFreeze - 1.9 _d 0
200 oceTfreez = celsius2K - 1.9 _d 0
201 DO J=1,NGP
202 sti1(J,myThid) = sst1(J,myThid)
203 IF ( oice1(J) .GT. 1. _d -2 ) THEN
204 sst1(J,myThid) = MAX(sst1(J,myThid),oceTfreez)
205 sti1(J,myThid) = sst1(J,myThid)
206 & +(sti1(J,myThid)-sst1(J,myThid))/oice1(J)
207 ELSE
208 oice1(J) = 0. _d 0
209 ENDIF
210 ENDDO
211 ELSE
212 DO J=1,NGP
213 sti1(J,myThid) = sst1(J,myThid)
214 ENDDO
215 ENDIF
216
217 C- Surface Albedo : (from F.M. FORDATE S/R)
218 c_FM DALB=ALBICE-ALBSEA
219 RSD=1. _d 0/SDALB
220 ALPHA= 2. _d 0*PI*(TYEAR+10. _d 0/365. _d 0)
221 RZEN = COS(ALPHA) * ( -23.45 _d 0 * deg2rad)
222 ZC = COS(RZEN)
223 ZS = SIN(RZEN)
224 DO j=1,sNy
225 DO i=1,sNx
226 c_FM SNOWC=MIN(1.,RSD*SNOW1(I,J))
227 c_FM ALBL=ALB0(I,J)+MAX(ALBSN-ALB0(I,J),0.0)*SNOWC
228 c_FM ALBS=ALBSEA+DALB*OICE1(I,J)
229 c_FM ALB1(I,J)=FMASK1(I,J)*ALBL+FMASK0(I,J)*ALBS
230 I2 = i+(j-1)*sNx
231 alb_land = aim_albedo(i,j,bi,bj)
232 & + MAX( 0. _d 0, ALBSN-aim_albedo(i,j,bi,bj) )
233 & *MIN( 1. _d 0, RSD*snow1(I2))
234 c alb_sea = ALBSEA + DALB*oice1(I2)
235 c alb1(I2,0,myThid) = alb_sea
236 c & + (alb_land - alb_sea)*fMask1(I2,1,myThid)
237 ALBSEA1 = ALBSEA
238 IF ( aim_selectOceAlbedo .EQ. 1) THEN
239 SJ = SIN(yC(i,j,bi,bj) * deg2rad)
240 CJ = COS(yC(i,j,bi,bj) * deg2rad)
241 TMPA = SJ*ZS
242 TMPB = CJ*ZC
243 TMPL = -TMPA/TMPB
244 IF (TMPL .GE. 1.0 _d 0) THEN
245 CZEN = 0.0 _d 0
246 ELSEIF (TMPL .LE. -1.0 _d 0) THEN
247 CZEN = (2.0 _d 0)*TMPA*PI
248 CZEN2= PI*((2.0 _d 0)*TMPA*TMPA + TMPB*TMPB)
249 CZEN = CZEN2/CZEN
250 ELSE
251 hlim = ACOS(TMPL)
252 CZEN = 2.0 _d 0*(TMPA*hlim + TMPB*SIN(hlim))
253 CZEN2= 2.0 _d 0*TMPA*TMPA*hlim
254 & + 4.0 _d 0*TMPA*TMPB*SIN(hlim)
255 & + TMPB*TMPB*( hlim + 0.5 _d 0*SIN(2.0 _d 0*hlim) )
256 CZEN = CZEN2/CZEN
257 ENDIF
258 ALBSEA1 = ( ( 2.6 _d 0 / (CZEN**(1.7 _d 0) + 0.065 _d 0) )
259 & + ( 15. _d 0 * (CZEN-0.1 _d 0) * (CZEN-0.5 _d 0)
260 & * (CZEN-1.0 _d 0) ) ) / 100.0 _d 0
261 ENDIF
262 alb1(I2,1,myThid) = alb_land
263 C_DE alb1(I2,2,myThid) = ALBSEA
264 alb1(I2,2,myThid) = 0.5 _d 0 * ALBSEA
265 & + 0.5 _d 0 * ALBSEA1
266 alb1(I2,3,myThid) = ALBICE
267 ENDDO
268 ENDDO
269
270 C-- else aim_useFMsurfBC
271 ELSE
272 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
273
274 C- safer to initialise output argument aim_sWght0,1
275 C even if they are not used when aim_useFMsurfBC=F
276 aim_sWght1 = 0. _d 0
277 aim_sWght0 = 1. _d 0
278
279 C- Set surface forcing fields needed by atmos. physics package
280 C 1. Albedo (between 0-1)
281 C 2. Sea Surface temperatures (in situ Temp. [K])
282 C 3. Land Surface temperatures (in situ Temp. [K])
283 C 4. Soil moisture (between 0-1)
284 C Snow depth, Sea Ice (<- no need for now)
285
286 C Set surface albedo data (in [0,1]) from aim_albedo to alb1 :
287 IF (aim_useMMsurfFc) THEN
288 DO j=1,sNy
289 DO i=1,sNx
290 I2 = i+(j-1)*sNx
291 alb1(I2,1,myThid) = aim_albedo(i,j,bi,bj)
292 alb1(I2,2,myThid) = aim_albedo(i,j,bi,bj)
293 alb1(I2,3,myThid) = aim_albedo(i,j,bi,bj)
294 ENDDO
295 ENDDO
296 ELSE
297 DO j=1,sNy
298 DO i=1,sNx
299 I2 = i+(j-1)*sNx
300 alb1(I2,1,myThid) = 0.
301 alb1(I2,2,myThid) = 0.
302 alb1(I2,3,myThid) = 0.
303 ENDDO
304 ENDDO
305 ENDIF
306 C Set surface temperature data from aim_S/LSurfTemp to sst1 & stl1 :
307 IF (aim_useMMsurfFc) THEN
308 DO j=1,sNy
309 DO i=1,sNx
310 I2 = i+(j-1)*sNx
311 sst1(I2,myThid) = aim_sst0(i,j,bi,bj)
312 stl1(I2,myThid) = aim_sst0(i,j,bi,bj)
313 sti1(I2,myThid) = aim_sst0(i,j,bi,bj)
314 ENDDO
315 ENDDO
316 ELSE
317 DO j=1,sNy
318 DO i=1,sNx
319 I2 = i+(j-1)*sNx
320 sst1(I2,myThid) = 300.
321 stl1(I2,myThid) = 300.
322 sti1(I2,myThid) = 300.
323 C_EqCh: start
324 sst1(I2,myThid) = 280.
325 & +20. _d 0 *exp( -((xC(i,j,bi,bj)-xBump)/dxBump)**2
326 & -((yC(i,j,bi,bj)-yBump)/dyBump)**2 )
327 stl1(I2,myThid) = sst1(I2,myThid)
328 sti1(I2,myThid) = sst1(I2,myThid)
329 C_EqCh: end
330 ENDDO
331 ENDDO
332 C_EqCh: start
333 IF (myIter.EQ.nIter0) THEN
334 WRITE(suff,'(I10.10)') myIter
335 CALL AIM_WRITE_PHYS( 'aim_SST.', suff, 1,sst1,
336 & 1, bi, bj, 1, myIter, myThid )
337 ENDIF
338 C_EqCh: end
339 ENDIF
340
341 C- Set soil water availability (in [0,1]) from aim_sw10 to soilw1 :
342 IF (aim_useMMsurfFc) THEN
343 DO j=1,sNy
344 DO i=1,sNx
345 I2 = i+(j-1)*sNx
346 soilw1(I2,myThid) = aim_sw10(i,j,bi,bj)
347 ENDDO
348 ENDDO
349 ELSE
350 DO j=1,sNy
351 DO i=1,sNx
352 I2 = i+(j-1)*sNx
353 soilw1(I2,myThid) = 0.
354 ENDDO
355 ENDDO
356 ENDIF
357
358 C- Set Snow depth and Sea Ice
359 C (not needed here since albedo is loaded from file)
360 DO j=1,sNy
361 DO i=1,sNx
362 I2 = i+(j-1)*sNx
363 oice1(I2) = 0.
364 snow1(I2) = 0.
365 ENDDO
366 ENDDO
367
368 C-- endif/else aim_useFMsurfBC
369 ENDIF
370
371 #ifdef COMPONENT_MODULE
372 IF ( useCoupler ) THEN
373 C-- take surface data from the ocean component
374 C to replace MxL fields (if use sea-ice) or directly AIM SST
375 CALL ATM_APPLY_IMPORT(
376 I aim_landFr,
377 U sst1(1,myThid), oice1,
378 I myTime, myIter, bi, bj, myThid )
379 ENDIF
380 #endif /* COMPONENT_MODULE */
381
382 #ifdef ALLOW_AIM_CO2
383 DO j=1,sNy
384 DO i=1,sNx
385 I2 = i+(j-1)*sNx
386 aim_CO2(I2,myThid)= atm_pCO2
387 ENDDO
388 ENDDO
389 #ifdef ALLOW_DIAGNOSTICS
390 IF ( useDiagnostics ) THEN
391 pCO2scl = 1. _d 6
392 CALL DIAGNOSTICS_SCALE_FILL( aim_CO2(1,myThid), pCO2scl, 1,
393 & 'aim_pCO2', 1, 1, 3, bi, bj, myThid )
394 ENDIF
395 #endif /* ALLOW_DIAGNOSTICS */
396 #endif /* ALLOW_AIM_CO2 */
397
398 #ifdef ALLOW_LAND
399 IF (useLand) THEN
400 C- Use land model output instead of prescribed Temp & moisture
401 CALL AIM_LAND2AIM(
402 I aim_landFr, aim_veget, aim_albedo, snow1,
403 U stl1(1,myThid), soilw1(1,myThid), alb1(1,1,myThid),
404 I myTime, myIter, bi, bj, myThid )
405 ENDIF
406 #endif /* ALLOW_LAND */
407
408 #ifdef ALLOW_THSICE
409 IF (useThSIce) THEN
410 C- Use thermo. sea-ice model output instead of prescribed Temp & albedo
411 CALL AIM_SICE2AIM(
412 I aim_landFr,
413 U sst1(1,myThid), oice1,
414 O sti1(1,myThid), alb1(1,3,myThid),
415 I myTime, myIter, bi, bj, myThid )
416 ENDIF
417 #endif /* ALLOW_THSICE */
418
419 C-- set the sea-ice & open ocean fraction :
420 DO J=1,NGP
421 fMask1(J,3,myThid) =(1. _d 0 - fMask1(J,1,myThid))
422 & *oice1(J)
423 fMask1(J,2,myThid) = 1. _d 0 - fMask1(J,1,myThid)
424 & - fMask1(J,3,myThid)
425 ENDDO
426
427 C-- set the mean albedo :
428 DO J=1,NGP
429 alb1(J,0,myThid) = fMask1(J,1,myThid)*alb1(J,1,myThid)
430 & + fMask1(J,2,myThid)*alb1(J,2,myThid)
431 & + fMask1(J,3,myThid)*alb1(J,3,myThid)
432 ENDDO
433
434 C-- initialize surf. temp. change to zero:
435 DO k=1,3
436 DO J=1,NGP
437 dTsurf(J,k,myThid) = 0.
438 ENDDO
439 ENDDO
440
441 IF (.NOT.aim_splitSIOsFx) THEN
442 DO J=1,NGP
443 fMask1(J,3,myThid) = 0. _d 0
444 fMask1(J,2,myThid) = 1. _d 0 - fMask1(J,1,myThid)
445 ENDDO
446 ENDIF
447
448 #endif /* ALLOW_AIM */
449
450 RETURN
451 END

  ViewVC Help
Powered by ViewVC 1.1.22