/[MITgcm]/MITgcm_contrib/darwin2/pkg/darwin/dic_surfforcing_init.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/darwin/dic_surfforcing_init.F

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


Revision 1.2 - (hide annotations) (download)
Wed Apr 20 19:19:27 2011 UTC (14 years, 8 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt62w_20110426
Changes since 1.1: +4 -1 lines
o add fix to two bugs in handling water vapour pressure in carbonate chemisty
  and air-sea fluxes. Found by Val Bennington and Galen McKinley

1 jahn 1.1 #include "CPP_OPTIONS.h"
2     #include "PTRACERS_OPTIONS.h"
3     #include "DARWIN_OPTIONS.h"
4    
5     #ifdef ALLOW_PTRACERS
6     #ifdef ALLOW_DARWIN
7    
8     #ifdef ALLOW_CARBON
9    
10     CBOP
11     C !ROUTINE: DIC_SURFFORCING_INIT
12    
13     C !INTERFACE: ==========================================================
14     SUBROUTINE DIC_SURFFORCING_INIT(
15     I myThid)
16    
17     C !DESCRIPTION:
18     C Calculate first guess of pH
19    
20     C !USES: ===============================================================
21     IMPLICIT NONE
22     #include "SIZE.h"
23     #include "DYNVARS.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "GRID.h"
27     #include "FFIELDS.h"
28     #include "PTRACERS_SIZE.h"
29     #include "PTRACERS_PARAMS.h"
30     #include "PTRACERS_FIELDS.h"
31     #include "DARWIN_SIZE.h"
32     #include "DARWIN_IO.h"
33     #include "DARWIN_FLUX.h"
34     #include "DIC_ATMOS.h"
35    
36     C !INPUT PARAMETERS: ===================================================
37     C myThid :: thread number
38     INTEGER myThid
39    
40    
41     C !LOCAL VARIABLES: ====================================================
42     INTEGER i,j, k, kLev, it
43     INTEGER intime0,intime1
44     _RL otime
45     _RL aWght,bWght,rdt
46     INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
47     C Number of iterations for pCO2 solvers...
48     C Solubility relation coefficients
49     C local variables for carbon chem
50     INTEGER iMin,iMax,jMin,jMax, bi, bj
51     _RL surfdic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52     _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53     _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54     _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55     _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56     INTEGER iprt,jprt
57     CHARACTER*(MAX_LEN_MBUF) msgBuf
58     INTEGER iUnit
59     CHARACTER*(MAX_LEN_FNAM) iName
60     integer ilo,ihi
61     integer ilnblnk,ifnblnk
62     external ilnblnk,ifnblnk
63     LOGICAL pH_isLoaded
64     CEOP
65    
66     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
67    
68     kLev=1
69    
70     CALL DIC_ATMOS(0, startTime, nIter0, myThid )
71    
72     _BEGIN_MASTER(myThid)
73    
74     C set up coefficients for DIC chemistry
75     C define Schmidt no. coefficients for CO2
76     sca1 = 2073.1 _d 0
77     sca2 = -125.62 _d 0
78     sca3 = 3.6276 _d 0
79     sca4 = -0.043219 _d 0
80     C define Schmidt no. coefficients for O2
81     C based on Keeling et al [GBC, 12, 141, (1998)]
82     sox1 = 1638.0 _d 0
83     sox2 = -81.83 _d 0
84     sox3 = 1.483 _d 0
85     sox4 = -0.008004 _d 0
86    
87     C coefficients for determining saturation O2
88     oA0= 2.00907 _d 0
89     oA1= 3.22014 _d 0
90     oA2= 4.05010 _d 0
91     oA3= 4.94457 _d 0
92     oA4= -2.56847 _d -1
93     oA5= 3.88767 _d 0
94     oB0= -6.24523 _d -3
95     oB1= -7.37614 _d -3
96     oB2= -1.03410 _d -2
97     oB3= -8.17083 _d -3
98     oC0= -4.88682 _d -7
99     C Set other constant/flag
100    
101     #ifndef USE_ATMOSCO2
102    
103     if (dic_int1.eq.2) then
104     call mdsfindunit( iUnit, mythid )
105     open(UNIT=iUnit,FILE='co2atmos.dat',STATUS='old')
106     do k=1,dic_int2
107     read(iUnit,*) co2atmos(k)
108     print*,'co2atmos',co2atmos(k)
109     enddo
110     close(iUnit)
111     endif
112    
113     if (dic_int1.eq.3) then
114     write(iName,'(A,I10.10)') 'dic_atmos.',nIter0
115     ilo = ifnblnk(iName)
116     ihi = ilnblnk(iName)
117     call mdsfindunit( iUnit, mythid )
118     open(UNIT=iUnit,FILE=iname(ilo:ihi),STATUS='old')
119     read(iUnit,*) total_atmos_carbon_ini,
120     & atpco2_ini
121     close(iUnit)
122     endif
123    
124     #endif
125     _END_MASTER(myThid)
126    
127     ccccccccccccccccccccccccccccccccccccccccc
128     IF ( periodicExternalForcing ) THEN
129    
130    
131     rdt = 1. _d 0 / deltaTclock
132     nForcingPeriods = NINT(externForcingCycle/externForcingPeriod)
133     cswd QQ change for placement of chem forcing (ie. after timestep)
134     Imytm = NINT(startTime*rdt)
135     Ifprd = NINT(externForcingPeriod*rdt)
136     Ifcyc = NINT(externForcingCycle*rdt)
137     Iftm = MOD( Imytm+Ifcyc-Ifprd/2, Ifcyc)
138    
139     intime0 = 1 + INT(Iftm/Ifprd)
140     intime1 = 1 + MOD(intime0,nForcingPeriods)
141     c aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
142     aWght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
143     bWght = FLOAT( Ifprd )
144     aWght = aWght / bWght
145     bWght = 1. _d 0 - aWght
146    
147     _BARRIER
148     _BEGIN_MASTER(myThid)
149    
150     _END_MASTER(myThid)
151    
152     #ifdef ALLOW_OFFLINE
153     IF ( useOffLine ) THEN
154     otime=nIter0*deltaTclock
155     CALL OFFLINE_FIELDS_LOAD( otime, nIter0, myThid )
156     ENDIF
157     #endif
158    
159     c end periodicExternalForcing
160     ENDIF
161    
162     C =================================================================
163    
164     jMin=1
165     jMax=sNy
166     iMin=1
167     iMax=sNx
168    
169     DO bj=myByLo(myThid),myByHi(myThid)
170     DO bi=myBxLo(myThid),myBxHi(myThid)
171     DO j=1-OLy,sNy+OLy
172     DO i=1-Olx,sNx+OLx
173     pH(i,j,bi,bj) = 8. _d 0
174     ENDDO
175     ENDDO
176     ENDDO
177     ENDDO
178    
179     DO bj=myByLo(myThid),myByHi(myThid)
180     DO bi=myBxLo(myThid),myBxHi(myThid)
181     DO j=1-OLy,sNy+OLy
182     DO i=1-OLx,sNx+OLx
183     ak0(i,j,bi,bj)=0. _d 0
184     ak1(i,j,bi,bj)=0. _d 0
185     ak2(i,j,bi,bj)=0. _d 0
186     akw(i,j,bi,bj)=0. _d 0
187     akb(i,j,bi,bj)=0. _d 0
188     akf(i,j,bi,bj)=0. _d 0
189     ak1p(i,j,bi,bj)=0. _d 0
190     ak2p(i,j,bi,bj)=0. _d 0
191     ak3p(i,j,bi,bj)=0. _d 0
192     aksi(i,j,bi,bj)=0. _d 0
193 stephd 1.2 fugf(i,j,bi,bj)=0. _d 0
194 jahn 1.1 ff(i,j,bi,bj)=0. _d 0
195     ft(i,j,bi,bj)=0. _d 0
196     st(i,j,bi,bj)=0. _d 0
197     bt(i,j,bi,bj)=0. _d 0
198     ENDDO
199     ENDDO
200     ENDDO
201     ENDDO
202    
203     pH_isLoaded = .FALSE.
204     IF ( nIter0.GT.PTRACERS_Iter0 ) THEN
205     C Read pH from a pickup file if needed
206     CALL DIC_READ_PICKUP(
207     O pH_isLoaded,
208     I nIter0, myThid )
209     ENDIF
210    
211     DO bj=myByLo(myThid),myByHi(myThid)
212     DO bi=myBxLo(myThid),myBxHi(myThid)
213    
214     C determine inorganic carbon chem coefficients
215     DO j=jMin,jMax
216     DO i=iMin,iMax
217    
218     c use surface layer values and convert to mol/m3
219     surfdic(i,j) = max(10. _d 0 ,
220     & Ptracer(i,j,1,bi,bj,iDIC))*1e-3
221     & * maskC(i,j,kLev,bi,bj)
222     surfalk(i,j) = max(10. _d 0 ,
223     & Ptracer(i,j,1,bi,bj,iALK))*1e-3
224     & * maskC(i,j,kLev,bi,bj)
225     surfphos(i,j) = max(1. _d -10,
226     & Ptracer(i,j,1,bi,bj,iPO4))*1e-3
227     & * maskC(i,j,kLev,bi,bj)
228     surfsi(i,j) = max(1. _d -8,
229     & Ptracer(i,j,1,bi,bj,iSi))*1e-3
230     & * maskC(i,j,kLev,bi,bj)
231     surfsalt(i,j) = max(4. _d 0, salt(i,j,kLev,bi,bj))
232     c
233     WIND(i,j,bi,bj) = 5. _d 0*maskC(i,j,1,bi,bj)
234     AtmosP(i,j,bi,bj) = 1. _d 0*maskC(i,j,1,bi,bj)
235     ENDDO
236     ENDDO
237    
238     CALL CARBON_COEFFS(
239     I theta,salt,
240     I bi,bj,iMin,iMax,jMin,jMax,myThid)
241    
242     C====================================================================
243    
244     IF ( .NOT.pH_isLoaded ) THEN
245     C set guess of pH for first step here
246    
247     print*,'QQ: pCO2 approximation method'
248     c first approximation
249     C$TAF LOOP = parallel
250     DO j=jMin,jMax
251     C$TAF LOOP = parallel
252     DO i=iMin,iMax
253     IF ( maskC(i,j,kLev,bi,bj) .NE. 0. _d 0) THEN
254     C$TAF init dic_surf = static, 10
255     DO it=1,10
256     C$TAF STORE pH(i,j,bi,bj), PTR_CO2(i,j,kLev) = dic_surf
257     C$TAF STORE surfalk(i,j), surfphos(i,j), surfsi(i,j) = dic_surf
258     CALL CALC_PCO2_APPROX(
259     I theta(i,j,kLev,bi,bj),surfsalt(i,j),
260     I surfdic(i,j), surfphos(i,j),
261     I surfsi(i,j),surfalk(i,j),
262     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
263     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
264     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
265 stephd 1.2 I aksi(i,j,bi,bj),akf(i,j,bi,bj),
266     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
267     I ff(i,j,bi,bj),
268 jahn 1.1 I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
269     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),
270     I myThid )
271     ENDDO
272     ENDIF
273     ENDDO
274     ENDDO
275     iprt = MIN(20,sNx)
276     jprt = MIN(20,sNy)
277     print*,'QQ first guess pH', pH(iprt,jprt,bi,bj),
278     & theta(iprt,jprt,1,bi,bj), salt(iprt,jprt,1,bi,bj),
279     & surfdic(iprt,jprt), surfphos(iprt,jprt),
280     & surfsi(iprt,jprt),surfalk(iprt,jprt)
281    
282     ENDIF
283    
284     C end bi,bj loops
285     ENDDO
286     ENDDO
287    
288     RETURN
289     END
290     #endif /*ALLOW_CARBON*/
291    
292     #endif /*DARWIN*/
293     #endif /*ALLOW_PTRACERS*/
294     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22