/[MITgcm]/MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F

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


Revision 1.1 - (hide annotations) (download)
Wed Apr 13 18:56:26 2011 UTC (14 years, 3 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt63c_20111011, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63d_20111107, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_baseline, ctrb_darwin2_ckpt63n_20120604, ctrb_darwin2_ckpt63k_20120317, ctrb_darwin2_ckpt62z_20110622
darwin2 initial checkin

1 jahn 1.1 C $Header: /u/gcmpack/MITgcm_contrib/darwin/pkg/darwin/darwin_forcing.F,v 1.16 2009/03/10 20:44:30 stephd Exp $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5     #include "PTRACERS_OPTIONS.h"
6     #include "DARWIN_OPTIONS.h"
7    
8     #ifdef ALLOW_PTRACERS
9     #ifdef ALLOW_DARWIN
10     #ifdef ALLOW_QUOTA
11    
12     c=============================================================
13     c subroutine quota_forcing
14     c step forward bio-chemical tracers in time
15     C==============================================================
16     SUBROUTINE QUOTA_FORCING(
17     U Ptr,
18     I bi,bj,imin,imax,jmin,jmax,
19     I myTime,myIter,myThid)
20     #include "SIZE.h"
21     #include "EEPARAMS.h"
22     #include "PARAMS.h"
23     #include "GRID.h"
24     #include "PTRACERS_SIZE.h"
25     #include "PTRACERS_PARAMS.h"
26     #include "GCHEM.h"
27     #include "QUOTA_SIZE.h"
28     #include "QUOTA.h"
29     #include "DARWIN_IO.h"
30     #include "DYNVARS.h"
31     #ifdef USE_QSW
32     #include "FFIELDS.h"
33     #endif
34    
35     C === Global variables ===
36     c tracers
37     _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)
38     INTEGER bi,bj,imin,imax,jmin,jmax
39     _RL myTime
40     INTEGER myIter
41     INTEGER myThid
42    
43     C============== Local variables ============================================
44     c biomodel tracer arrays
45     _RL nutrient(iimax)
46     _RL biomass(iomax,npmax)
47     _RL orgmat(iomax-iChl,komax)
48     #ifdef FQUOTA
49     c iron partitioning
50     _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51     _RL freefu
52     _RL inputFel
53     #endif
54     c upstream arrays for sinking/swimming
55     _RL bioabove(iomax,npmax)
56     _RL biobelow(iomax,npmax)
57     _RL orgabove(iomax-iChl,komax)
58     c some working variables
59     _RL sumpy
60     _RL sumpyup
61     c light variables
62     _RL PAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63     _RL sfac(1-OLy:sNy+OLy)
64     _RL atten,lite
65     _RL newtime ! for sub-timestepping
66     _RL runtim ! time from tracer initialization
67    
68     #ifdef DAR_DIAG_DIVER
69     _RL Diver1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
70     _RL Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
71     _RL Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72     _RL Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73     #endif
74     c
75    
76     c some local variables
77     _RL Tlocal
78     _RL Slocal
79     _RL PARlocal
80     _RL dzlocal
81     _RL dtplankton
82     _RL PP
83     c local tendencies
84     _RL dbiomass(iomax,npmax)
85     _RL dorgmat(iomax-iChl,komax)
86     _RL dnutrient(iimax)
87     _RL tmp
88    
89     INTEGER bottom
90     INTEGER surface
91     INTEGER i,j,k,it, ktmp
92     INTEGER ii,io,jp,ko, jp2, jpsave
93     INTEGER place
94     INTEGER debug
95     CHARACTER*8 diagname
96    
97     c
98     c--------------------------------------------------
99     c initialise vatriables
100     DO j=1-OLy,sNy+OLy
101     DO i=1-OLx,sNx+OLx
102     do k=1,Nr
103     #ifdef FQUOTA
104     freefe(i,j,k) = 0.0 _d 0
105     # endif
106     PAR(i,j,k) = 0.0 _d 0
107     #ifdef DAR_DIAG_DIVER
108     Diver1(i,j,k) = 0.0 _d 0
109     Diver2(i,j,k) = 0.0 _d 0
110     Diver3(i,j,k) = 0.0 _d 0
111     Diver4(i,j,k) = 0.0 _d 0
112     #endif
113     c
114     enddo !k
115     ENDDO !i
116     ENDDO !j
117     c
118     c bio-chemical time loop
119     c--------------------------------------------------
120     DO it=1,nsubtime
121     c -------------------------------------------------
122     COJ cannot use dfloat because of adjoint
123     COJ division will be double precision anyway because of dTtracerLev
124     newtime=myTime-dTtracerLev(1)+
125     & float(it)*dTtracerLev(1)/float(nsubtime)
126     c print*,'it ',it,newtime,nsubtime,myTime
127     runtim=myTime-float(PTRACERS_Iter0)*dTtracerLev(1)
128    
129     #ifdef FQUOTA
130     c determine iron partitioning - solve for free iron
131     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
132     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
133     & myIter, mythid)
134     #endif
135    
136     c find light in each grid cell
137     c ---------------------------
138     c determine incident light
139     #ifndef READ_PAR
140     #ifdef USE_QSW
141     DO j=1-OLy,sNy+OLy
142     DO i=1-OLx,sNx+OLx
143     sur_par(i,j,bi,bj)=-parfrac*Qsw(i,j,bi,bj)*
144     & parconv*maskC(i,j,1,bi,bj)
145     ENDDO
146     ENDDO
147     #else
148     DO j=1-OLy,sNy+OLy
149     sfac(j)=0. _d 0
150     ENDDO
151     call darwin_insol(newTime,sfac,bj)
152     DO j=1-OLy,sNy+OLy
153     DO i=1-OLx,sNx+OLx
154     sur_par(i,j,bi,bj)=sfac(j)*maskC(i,j,1,bi,bj)/86400. _d 6
155     c if (i.eq.1.and.j.ge.1.and.j.le.sNy)
156     c & write(24,*) sur_par(i,j,bi,bj)
157     ENDDO
158     ENDDO
159     #endif
160     #endif
161    
162     C.................................................................
163     C.................................................................
164    
165    
166     DO j=1,sNy
167     DO i=1,sNx
168     c surface PAR
169     c take ice coverage into account
170     #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
171     COJ ice coverage already taken into account by seaice package
172     lite=sur_par(i,j,bi,bj)
173     #else
174     #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
175     c if using Qsw and seaice, then ice fraction is already
176     c taken into account
177     lite=sur_par(i,j,bi,bj)
178     #else
179     lite=sur_par(i,j,bi,bj)*(1. _d 0-fice(i,j,bi,bj))
180     #endif
181     #endif
182     atten = 0. _d 0
183     sumpy = 0. _d 0
184     c
185     c FOR EACH LAYER ...
186     do k= 1, NR
187     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
188     c ---------------------------------------------------------------------
189     c benw
190     c
191     c Fetch biomodel variables from ptr (ptracers)
192     c (making sure they are .ge. 0 - brute force)
193     c
194     c (set biomodel tendencies to zero, at the same time)
195     c
196     c *********************************************************************
197     place = 0
198     c Inorganic Nutrients
199     do ii=1,iimax
200     place = place + 1
201     c ambient nutrients for each element (1 to iimax)
202     nutrient(ii) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
203     dnutrient(ii) = 0. _d 0
204     enddo ! ii
205     c *********************************************************************
206     c Unicellular biomass (including chlorophyll biomass)
207     do io=1,iomax
208     do jp=1,npmax
209     place = place + 1
210     biomass(io,jp) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
211     ! biomasses above current layer for sinking
212     if (k.eq.1) then
213     bioabove(io,jp)=0. _d 0
214     endif
215     ! biomasses below current layer for swimming
216     if (k.eq.Nr) then
217     biobelow(io,jp)=0. _d 0
218     elseif (hFacC(i,j,k+1,bi,bj).eq.0. _d 0) then
219     biobelow(io,jp)=0. _d 0
220     else
221     biobelow(io,jp)=max(Ptr(i,j,k+1,bi,bj,place),0. _d 0)
222     endif
223     ! initialise biomass rate of change
224     dbiomass(io,jp) = 0. _d 0
225     enddo ! jp
226     enddo
227     c *********************************************************************
228     c Organic matter
229     do io=1,iomax-iChl
230     do ko=1,komax
231     c mass of element x for all OM classes
232     place = place + 1
233     orgmat(io,ko) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
234     ! biomasses above current layer for sinking
235     if (k.eq.1) then
236     orgabove(io,ko) = 0. _d 0
237     endif
238     #ifdef SQUOTA
239     if (ko.and.1.and.io.eq.iSili) then
240     place = place - 1
241     orgmat(iSili,1) = 0. _d 0
242     orgabove(iSili,1) = 0. _d 0
243     endif
244     #endif
245     dorgmat(io,ko) = 0. _d 0
246     enddo ! ko
247     enddo ! io
248     c *********************************************************************
249     c
250     c ---------------------------------------------------------------------
251    
252    
253     c find local light for level k
254     sumpyup = sumpy
255     sumpy = 0. _d 0
256     do jp=1,npmax
257     #ifndef GEIDER
258     ! sum nitrogen biomass
259     sumpy = sumpy + biomass(iNitr,jp)
260     #else
261     ! sum chlorophyll
262     sumpy = sumpy + biomass(iChlo,jp)
263     #endif
264     enddo
265    
266     atten= atten + (k_w + k_chl*sumpy)*5. _d -1*drF(k)
267     if (k.gt.1)then
268     atten = atten + (k_w+k_chl*sumpyup)*5. _d -1*drF(k-1)
269     endif
270     PAR(i,j,k) = lite*exp(-atten)
271     c
272     c Physical variables
273     PARlocal = PAR(i,j,k)
274     Tlocal = theta(i,j,k,bi,bj)
275     Slocal = salt(i,j,k,bi,bj)
276     c Free Iron
277     #ifdef FQUOTA
278     freefu = max(freefe(i,j,k),0. _d 0)
279     if (k.eq.1) then
280     inputFel = inputFe(i,j,bi,bj)
281     else
282     inputFel = 0. _d 0
283     endif
284     #endif
285     c Layer thickness
286     dzlocal = drF(k)*HFacC(i,j,k,bi,bj)
287     c
288     c set bottom=1.0 if the layer below is not ocean
289     ktmp=min(nR,k+1)
290     if(hFacC(i,j,ktmp,bi,bj).eq.0. _d 0.or.k.eq.Nr) then
291     bottom = 1
292     else
293     bottom = 0
294     endif
295     if (k.eq.1) then
296     surface = 1
297     else
298     surface = 0
299     endif
300    
301     c set other arguments to zero
302     debug=0
303    
304     if (debug.eq.7) print*,'Inorganic nutrients',nutrient
305     if (debug.eq.7) print*,'Plankton biomass', biomass
306     if (debug.eq.7) print*,'Organic nutrients',orgmat
307     if (debug.eq.8) print*,'k, PARlocal, dzlocal',
308     & k,PARlocal,dzlocal
309     c ---------------------------------------------------------------------
310     CALL QUOTA_PLANKTON(
311     I biomass, orgmat, nutrient,
312     O PP,
313     I bioabove, biobelow,
314     I orgabove,
315     #ifdef FQUOTA
316     I freefu, inputFel,
317     #endif
318     I PARlocal, Tlocal, Slocal,
319     I bottom, surface, dzlocal,
320     O dbiomass, dorgmat, dnutrient,
321     I debug,
322     I runtim,
323     I MyThid)
324     c ---------------------------------------------------------------------
325     c save un-updated biomass as layer above
326     do io=1,iomax
327     do jp=1,npmax
328     bioabove(io,jp)=biomass(io,jp)
329     enddo
330     if (io.ne.iChlo) then
331     do ko=1,komax
332     orgabove(io,ko)=orgmat(io,ko)
333     enddo
334     endif
335     enddo
336     c ---------------------------------------------------------------------
337     c now update main tracer arrays
338     c for timestep dtplankton
339     dtplankton = dTtracerLev(k)/float(nsubtime)
340     cccccccccccccccccccccccccccccccccccccccccccccccccccc
341     place = 0
342     cccccccccccccccccccccccccccccccccccccccccccccccccccc
343     c Inorganic nutrients
344     do ii=1,iimax
345     place = place + 1
346     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
347     & + dtplankton*dnutrient(ii)
348     enddo ! ii
349     cccccccccccccccccccccccccccccccccccccccccccccccccccc
350     c Biomass
351     do io=1,iomax
352     do jp=1,npmax
353     place = place + 1
354     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
355     & + dtplankton*dbiomass(io,jp)
356     if (pft(jp).eq.6.and.io.eq.iChlo) then
357     Ptr(i,j,k,bi,bj,place) = 0. _d 0
358     endif
359     enddo ! jp
360     enddo ! io
361     ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
362     c Organic matter
363     do io=1,iomax-iChl
364     do ko=1,komax
365     if (ko.ne.1.or.io.ne.iSili) then
366     place = place + 1
367     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
368     & + dtplankton*dorgmat(io,ko)
369     endif
370     enddo ! ko
371     enddo ! io
372     ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
373     c
374     PPave(i,j,k,bi,bj) = PPave(i,j,k,bi,bj)+
375     & PP*dtplankton
376     PARave(i,j,k,bi,bj) = PARave(i,j,k,bi,bj)+
377     & PARlocal * dtplankton
378     c
379     #ifdef ALLOW_TIMEAVE
380     c
381     #ifdef DAR_DIAG_DIVER
382     Diver1ave(i,j,k,bi,bj)=Diver1ave(i,j,k,bi,bj)+
383     & Diver1(i,j,k)*dtplankton
384     Diver2ave(i,j,k,bi,bj)=Diver2ave(i,j,k,bi,bj)+
385     & Diver2(i,j,k)*dtplankton
386     Diver3ave(i,j,k,bi,bj)=Diver3ave(i,j,k,bi,bj)+
387     & Diver3(i,j,k)*dtplankton
388     Diver4ave(i,j,k,bi,bj)=Diver4ave(i,j,k,bi,bj)+
389     & Diver4(i,j,k)*dtplankton
390     #endif
391     #endif
392     endif
393     c end if hFac>0
394     enddo ! k
395     c end layer loop
396     c
397     ENDDO ! i
398     ENDDO ! j
399     c
400     c
401     COJ fill diagnostics
402     #ifdef ALLOW_DIAGNOSTICS
403     IF ( useDiagnostics ) THEN
404     diagname = ' '
405     do jp=1,npmax
406     WRITE(diagname,'(A8)') 'dCHL',jp,' '
407     CALL DIAGNOSTICS_FILL
408     & (dCHLarr(1-Olx,1-Oly,1,jp),diagname,0,Nr,2,bi,bj,myThid)
409     do ii=1,iimax
410     WRITE(diagname,'(A8)') 'PP',ii,jp,' '
411     CALL DIAGNOSTICS_FILL
412     & (PParr(1-Olx,1-Oly,1,ii,jp),diagname,0,Nr,2,bi,bj,myThid)
413     enddo
414     enddo
415     c
416     WRITE(diagname,'(A8)') 'PAR '
417     CALL DIAGNOSTICS_FILL( PAR(1-Olx,1-Oly,1), diagname,
418     & 0,Nr,2,bi,bj,myThid )
419     #ifdef DAR_DIAG_DIVER
420     WRITE(diagname,'(A8)') 'Diver1 '
421     CALL DIAGNOSTICS_FILL( Diver1(1-Olx,1-Oly,1), diagname,
422     & 0,Nr,2,bi,bj,myThid )
423     WRITE(diagname,'(A8)') 'Diver2 '
424     CALL DIAGNOSTICS_FILL( Diver2(1-Olx,1-Oly,1), diagname,
425     & 0,Nr,2,bi,bj,myThid )
426     WRITE(diagname,'(A8)') 'Diver3 '
427     CALL DIAGNOSTICS_FILL( Diver3(1-Olx,1-Oly,1), diagname,
428     & 0,Nr,2,bi,bj,myThid )
429     WRITE(diagname,'(A8)') 'Diver4 '
430     CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
431     & 0,Nr,2,bi,bj,myThid )
432     #endif
433     ENDIF
434     #endif
435     COJ
436    
437     #ifdef FQUOTA
438     c determine iron partitioning - solve for free iron
439     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
440     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
441     & myIter, mythid)
442     #endif
443    
444     c
445     #ifdef ALLOW_TIMEAVE
446     c save averages
447     do k=1,nR
448     dar_timeave(bi,bj,k) = dar_timeave(bi,bj,k)
449     & + dtplankton
450     enddo
451     #endif
452     c
453     c -----------------------------------------------------
454     ENDDO ! it
455     c -----------------------------------------------------
456     c end of bio-chemical time loop
457     c
458     RETURN
459     END
460    
461     #endif /*ALLOW_QUOTA*/
462     #endif /*ALLOW_DARWIN*/
463     #endif /*ALLOW_PTRACERS*/
464    
465     C============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22