--- MITgcm/verification/global_with_CFC11/code1x1/external_forcing_tr.F 2005/08/25 16:22:17 1.1 +++ MITgcm/verification/global_with_CFC11/code1x1/external_forcing_tr.F 2005/08/25 18:21:17 1.1.2.2 @@ -0,0 +1,295 @@ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/global_with_CFC11/code1x1/Attic/external_forcing_tr.F,v 1.1.2.2 2005/08/25 18:21:17 dimitri Exp $ +C $Name: $ + +#include "CPP_OPTIONS.h" + +CBOP +C !ROUTINE: EXTERNAL_FORCING_TR +C !INTERFACE: + SUBROUTINE EXTERNAL_FORCING_TR( + I iMin, iMax, jMin, jMax,bi,bj,kLev, + I myTime,myThid) + +C !DESCRIPTION: \bv +C *==========================================================* +C | S/R EXTERNAL_FORCING_TR +C | o Contains problem specific forcing for tracer Tr1. +C *==========================================================* +C | Adds terms to gTr1 for forcing by external sources. +C | This routine is hardcoded for OCMIP CFC-11 experiment. +C | It assumes that myTime=0 on January 1 for FICE, XKW, +C | and PATM, and data.cal provides the date for cfc1112.atm. +C | See the OCMIP-2 CFC HOWTO for details. +C *==========================================================* +C \ev + +C !USES: + IMPLICIT NONE +C == Global data == +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "DYNVARS.h" +#include "TR1.h" + + _RL sc_cfc, sol_cfc + EXTERNAL sc_cfc, sol_cfc + +C !INPUT/OUTPUT PARAMETERS: +C == Routine arguments == +C iMin - Working range of tile for applying forcing. +C iMax +C jMin +C jMax +C kLev + INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj + _RL myTime + INTEGER myThid + +C !LOCAL VARIABLES: +C == Local variables == +C i,j :: Loop counters +C aWght, bWght :: Interpolation weights + INTEGER i,j,intime0,intime1,nForcingPeriods + INTEGER currentdate(4), tempdate(4), timediff(4) + _RL aWght,bWght,rdt,timeint + _RL fprd,fcyc,ftm,ftmold + _RL KW,CSAT,ALPHA,PCFC,PCFC1,PCFC2,TMP1,TMP2 +CEOP + + IF ( kLev .EQ. 1 ) THEN + +C Find records and compute interpolation weights + fprd = 2629800 + fcyc = 31557600 + nForcingPeriods = fcyc / fprd + ftm = mod( myTime + fcyc - fprd / 2 , fcyc ) + intime0 = int( ftm / fprd ) + intime1 = mod( intime0 + 1, nForcingPeriods ) + aWght = ( ftm - fprd * intime0 ) / ( fprd ) + bWght = 1. - aWght + intime0 = intime0 + 1 + intime1 = intime1 + 1 + +C Now calculate whether it is time to update the forcing arrays + ftmold = mod( myTime - deltaTclock + fcyc - fprd / 2 , fcyc ) + IF ( + & (int(ftmold/fprd)+1) .NE. intime0 + & .OR. myTime .EQ. startTime + & ) THEN +C If the above condition is met then we need to read in +C data for the period ahead and the period behind myTime. + _BEGIN_MASTER(myThid) + WRITE(*,*) + & 'S/R EXTERNAL_FORCING_TR: Reading new data',myTime + +#define CFC_USE_INTERP +C-- CFC_USE_INTERP option is useful for high-resolution integrations. +C For low-resolution, less that 1-deg, it's best to generate files +C separately because of sea-ice fraction. +C Caution: CFC_USE_INTERP as used is not thread-safe. +#ifdef CFC_USE_INTERP + CALL EXF_INTERP( FiceFile,readBinaryPrec,fice0,intime0, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) + CALL EXF_INTERP( FiceFile,readBinaryPrec,fice1,intime1, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) + CALL EXF_INTERP( XkwFile,readBinaryPrec,xkw0,intime0, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) + CALL EXF_INTERP( XkwFile,readBinaryPrec,xkw1,intime1, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) + CALL EXF_INTERP( PatmFile,readBinaryPrec,patm0,intime0, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) + CALL EXF_INTERP( PatmFile,readBinaryPrec,patm1,intime1, + & xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid ) +#else + CALL MDSREADFIELD ( FiceFile, readBinaryPrec, + & 'RS', 1, fice0, intime0, myThid ) + CALL MDSREADFIELD ( FiceFile, readBinaryPrec, + & 'RS', 1, fice1, intime1, myThid ) + CALL MDSREADFIELD ( XkwFile, readBinaryPrec, + & 'RS', 1, xkw0, intime0, myThid ) + CALL MDSREADFIELD ( XkwFile, readBinaryPrec, + & 'RS', 1, xkw1, intime1, myThid ) + CALL MDSREADFIELD ( PatmFile, readBinaryPrec, + & 'RS', 1, patm0, intime0, myThid ) + CALL MDSREADFIELD ( PatmFile, readBinaryPrec, + & 'RS', 1, patm1, intime1, myThid ) +#endif + + _END_MASTER(myThid) + _EXCH_XY_R4(fice0 , myThid ) + _EXCH_XY_R4(fice1 , myThid ) + _EXCH_XY_R4(xkw0 , myThid ) + _EXCH_XY_R4(xkw1 , myThid ) + _EXCH_XY_R4(patm0 , myThid ) + _EXCH_XY_R4(patm1 , myThid ) + ENDIF + +C-- Interpolate in time + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + fice(i,j,bi,bj) = bWght*fice0(i,j,bi,bj) + + & aWght*fice1(i,j,bi,bj) + if( fice(i,j,bi,bj) .LT. 0.2 ) fice(i,j,bi,bj) = 0.0 + xkw(i,j,bi,bj) = bWght* xkw0(i,j,bi,bj) + + & aWght* xkw1(i,j,bi,bj) + patm(i,j,bi,bj) = bWght*patm0(i,j,bi,bj) + + & aWght*patm1(i,j,bi,bj) + ENDDO + ENDDO + +C Find records and interpolation weights for PCFC + call cal_FullDate( 19310101, 0, tempdate, mythid ) + call cal_GetDate( 0, mytime, currentdate, mythid ) + call cal_SubDates( currentdate, tempdate, timediff, mythid ) + if( timediff(1) .GT. 0 ) then + call cal_ToSeconds( timediff, timeint, mythid ) + else + timeint=0 + endif + fprd = 31557600 + ftm = timeint + fprd / 2 + intime0 = int( ftm / fprd ) + intime1 = intime0 + 1 + aWght = ( ftm - fprd * intime0 ) / ( fprd ) + bWght = 1. - aWght + intime0 = intime0 + 1930 + intime1 = intime1 + 1930 + PCFC1 = bWght * CFCp11(intime0,1) + aWght * CFCp11(intime1,1) + PCFC2 = bWght * CFCp11(intime0,2) + aWght * CFCp11(intime1,2) + +C-- Forcing term: add tracer in top-layer +C TR1 is tracer concentration in mol/m^3 +C OCMIP formula provides air-sea flux F in mol/m^2/s +C GTR1 = F / drF(1) in mol/m^3/s +C + DO j=jMin,jMax + DO i=iMin,iMax + TMP1 = theta(i,j,1,bi,bj) + TMP2 = salt(i,j,1,bi,bj) + KW = (1.-fice(i,j,bi,bj)) * xkw(i,j,bi,bj) * + & (sc_cfc(TMP1,11)/660)**(-1/2) + PCFC = pCFCw1(I,J,bi,bj) * PCFC1 + + & pCFCw2(I,J,bi,bj) * PCFC2 + ALPHA = sol_cfc(TMP1,TMP2,11) + CSAT = ALPHA * PCFC * patm(i,j,bi,bj) + TMP1 = KW * ( CSAT - Tr1(i,j,1,bi,bj) ) + surfaceTendencyTr1(i,j,bi,bj) = + & TMP1 * recip_drF(1) * recip_hFacC(i,j,1,bi,bj) + ENDDO + ENDDO + + ENDIF + + RETURN + END + +C==================================================================== + + _RL FUNCTION sc_cfc(t,kn) + +c--------------------------------------------------- +c CFC 11 and 12 schmidt number +c as a fonction of temperature. +c +c ref: Zheng et al (1998), JGR, vol 103,No C1 +c +c t: temperature (degree Celcius) +c kn: = 11 for CFC-11, 12 for CFC-12 +c +c J-C Dutay - LSCE +c--------------------------------------------------- + + IMPLICIT NONE + INTEGER kn + _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12) + _RL t +c +c coefficients with t in degre Celcius +c ------------------------------------ + a1(11) = 3501.8 + a2(11) = -210.31 + a3(11) = 6.1851 + a4(11) = -0.07513 +c + a1(12) = 3845.4 + a2(12) = -228.95 + a3(12) = 6.1908 + a4(12) = -0.067430 +c + + sc_cfc = a1(kn) + a2(kn) * t + a3(kn) *t*t + & + a4(kn) *t*t*t + + RETURN + END + +C==================================================================== + + _RL FUNCTION sol_cfc(pt,ps,kn) + +c------------------------------------------------------------------- +c +c CFC 11 and 12 Solubilities in seawater +c ref: Warner & Weiss (1985) , Deep Sea Research, vol32 +c +c pt: temperature (degre Celcius) +c ps: salinity (o/oo) +c kn: 11 = CFC-11, 12 = CFC-12 +c sol_cfc: in mol/m3/pptv +c 1 pptv = 1 part per trillion = 10^-12 atm = 1 picoatm +c +c J-C Dutay - LSCE +c------------------------------------------------------------------- + + _RL pt, ps,ta,d + _RL a1 ( 11: 12), a2 ( 11: 12), a3 ( 11: 12), a4 ( 11: 12) + _RL b1 ( 11: 12), b2 ( 11: 12), b3 ( 11: 12) + + INTEGER kn + +cc +cc coefficient for solubility in mol/l/atm +cc ---------------------------------------- +c +c for CFC 11 +c ---------- + a1 ( 11) = -229.9261 + a2 ( 11) = 319.6552 + a3 ( 11) = 119.4471 + a4 ( 11) = -1.39165 + b1 ( 11) = -0.142382 + b2 ( 11) = 0.091459 + b3 ( 11) = -0.0157274 +c +c for CFC/12 +c ---------- + a1 ( 12) = -218.0971 + a2 ( 12) = 298.9702 + a3 ( 12) = 113.8049 + a4 ( 12) = -1.39165 + b1 ( 12) = -0.143566 + b2 ( 12) = 0.091015 + b3 ( 12) = -0.0153924 +c + + ta = ( pt + 273.16)* 0.01 + d = ( b3 ( kn)* ta + b2 ( kn))* ta + b1 ( kn) +c +c + sol_cfc + $ = exp ( a1 ( kn) + $ + a2 ( kn)/ ta + $ + a3 ( kn)* log ( ta ) + $ + a4 ( kn)* ta * ta + ps* d ) +c +c conversion from mol/(l * atm) to mol/(m^3 * atm) +c ------------------------------------------------ + sol_cfc = 1000. * sol_cfc +c +c conversion from mol/(m^3 * atm) to mol/(m3 * pptv) +c -------------------------------------------------- + sol_cfc = 1.0e-12 * sol_cfc + + END