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

Diff of /MITgcm_contrib/darwin2/pkg/darwin/dic_surfforcing.F

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

revision 1.2 by stephd, Wed Apr 20 19:19:27 2011 UTC revision 1.3 by stephd, Mon May 9 20:17:39 2011 UTC
# Line 31  C !USES: =============================== Line 31  C !USES: ===============================
31  #include "DARWIN_SIZE.h"  #include "DARWIN_SIZE.h"
32  #include "DARWIN_IO.h"  #include "DARWIN_IO.h"
33  #include "DARWIN_FLUX.h"  #include "DARWIN_FLUX.h"
34    #ifdef USE_EXFWIND
35    #include "EXF_FIELDS.h"
36    #endif
37    
38  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
39  C  myThid               :: thread number  C  myThid               :: thread number
# Line 64  C local variables for carbon chem Line 67  C local variables for carbon chem
67        _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68        _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69        _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70          _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71  #ifdef ALLOW_OLD_VIRTUALFLUX  #ifdef ALLOW_OLD_VIRTUALFLUX
72        _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73  #endif  #endif
# Line 87  C ====================================== Line 91  C ======================================
91  C determine inorganic carbon chem coefficients  C determine inorganic carbon chem coefficients
92          DO j=jmin,jmax          DO j=jmin,jmax
93           DO i=imin,imax           DO i=imin,imax
94               surfdic(i,j) = max(10. _d 0 , PTR_DIC(i,j))*1e-3  c put bounds on tracers so pH solver doesn't blow up
95                 surfdic(i,j) =
96         &           max(400. _d 0 , min(4000. _d 0, PTR_DIC(i,j)))*1e-3
97       &                          * maskC(i,j,kLev,bi,bj)       &                          * maskC(i,j,kLev,bi,bj)
98               surfalk(i,j) = max(10. _d 0 , PTR_ALK(i,j))*1e-3               surfalk(i,j) =
99         &           max(400. _d 0 , min(4000. _d 0, PTR_ALK(i,j)))*1e-3
100       &                          * maskC(i,j,kLev,bi,bj)       &                          * maskC(i,j,kLev,bi,bj)
101               surfphos(i,j)  = max(1. _d -10, PTR_PO4(i,j))*1e-3               surfphos(i,j)  =
102         &           max(1. _d -10, min(10. _d 0, PTR_PO4(i,j)))*1e-3
103       &                          * maskC(i,j,kLev,bi,bj)       &                          * maskC(i,j,kLev,bi,bj)
104               surfsi(i,j)   = max(1. _d -8, PTR_SIL(i,j))*1e-3               surfsi(i,j)   =
105         &           max(1. _d -8, min(500. _d 0, PTR_SIL(i,j)))*1e-3
106       &                          * maskC(i,j,kLev,bi,bj)       &                          * maskC(i,j,kLev,bi,bj)
107               surfsalt(i,j) = max(4. _d 0, salt(i,j,kLev,bi,bj))               surfsalt(i,j) =
108         &           max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
109                 surftemp(i,j) =
110         &            max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
111           ENDDO           ENDDO
112          ENDDO          ENDDO
113    
114           CALL CARBON_COEFFS(           CALL CARBON_COEFFS(
115       I                       theta,salt,       I                       surftemp,surfsalt,
116       I                       bi,bj,iMin,iMax,jMin,jMax,myThid)       I                       bi,bj,iMin,iMax,jMin,jMax,myThid)
117  C====================================================================  C====================================================================
118    
# Line 119  C       average AtmosP is about 1 Atm. Line 131  C       average AtmosP is about 1 Atm.
131    
132  C Pre-compute part of exchange coefficient: pisvel*(1-fice)  C Pre-compute part of exchange coefficient: pisvel*(1-fice)
133  C Schmidt number is accounted for later  C Schmidt number is accounted for later
134    #ifdef USE_EXFWIND
135                  pisvel(i,j)=0.337 _d 0 *wspeed(i,j,bi,bj)**2/3.6 _d 5
136    #else
137                pisvel(i,j)=0.337 _d 0 *wind(i,j,bi,bj)**2/3.6 _d 5                pisvel(i,j)=0.337 _d 0 *wind(i,j,bi,bj)**2/3.6 _d 5
138    #endif
139                Kwexch_Pre(i,j,bi,bj) = pisvel(i,j)                Kwexch_Pre(i,j,bi,bj) = pisvel(i,j)
140       &                              * (1. _d 0 - FIce(i,j,bi,bj))       &                              * (1. _d 0 - FIce(i,j,bi,bj))
141    
# Line 134  C$TAF LOOP = parallel Line 150  C$TAF LOOP = parallel
150    
151            IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN            IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
152              CALL CALC_PCO2_APPROX(              CALL CALC_PCO2_APPROX(
153       I        theta(i,j,kLev,bi,bj),surfsalt(i,j),       I        surftemp(i,j),surfsalt(i,j),
154       I        surfdic(i,j), surfphos(i,j),       I        surfdic(i,j), surfphos(i,j),
155       I        surfsi(i,j),surfalk(i,j),       I        surfsi(i,j),surfalk(i,j),
156       I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),       I        ak1(i,j,bi,bj),ak2(i,j,bi,bj),
# Line 160  C$TAF LOOP = parallel Line 176  C$TAF LOOP = parallel
176  C calculate SCHMIDT NO. for CO2  C calculate SCHMIDT NO. for CO2
177                SchmidtNoDIC(i,j) =                SchmidtNoDIC(i,j) =
178       &            sca1       &            sca1
179       &          + sca2 * theta(i,j,kLev,bi,bj)       &          + sca2 * surftemp(i,j)
180       &          + sca3 * theta(i,j,kLev,bi,bj)*theta(i,j,kLev,bi,bj)         &          + sca3 * surftemp(i,j)*surftemp(i,j)
181       &          + sca4 * theta(i,j,kLev,bi,bj)*theta(i,j,kLev,bi,bj)       &          + sca4 * surftemp(i,j)*surftemp(i,j)
182       &                *theta(i,j,kLev,bi,bj)       &                *surftemp(i,j)
183    c put positive bound on SCHMIT number (will go negative for temp>40)
184                  SchmidtNoDIC(i,j) = max(1. _d -2, SchmidtNoDIC(i,j))
185    
186  C Determine surface flux (FDIC)  C Determine surface flux (FDIC)
187  C first correct pCO2at for surface atmos pressure  C first correct pCO2at for surface atmos pressure
# Line 171  C first correct pCO2at for surface atmos Line 189  C first correct pCO2at for surface atmos
189       &          AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)       &          AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)
190    
191  C then account for Schmidt number  C then account for Schmidt number
192                Kwexch(i,j) = Kwexch_Pre(i,j,bi,bj)                  Kwexch(i,j) = Kwexch_Pre(i,j,bi,bj)
193       &                    / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)       &                    / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
194    
   
195  #ifdef WATERVAP_BUG  #ifdef WATERVAP_BUG
196  C Calculate flux in terms of DIC units using K0, solubility  C Calculate flux in terms of DIC units using K0, solubility
197  C Flux = Vp * ([CO2sat] - [CO2])  C Flux = Vp * ([CO2sat] - [CO2])
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22