pkg/bling/bling_light.F

c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
a284455135 Matt*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
c0d1c06c15 Matt*0005
                0006 CBOP
a284455135 Matt*0007       SUBROUTINE BLING_LIGHT(
e0f9a7ba0b Matt*0008      I           mld,
                0009      U           irr_inst, irr_eff,
                0010      I           bi, bj, imin, imax, jmin, jmax,
                0011      I           myTime, myIter, myThid)
                0012
c0d1c06c15 Matt*0013 C     =================================================================
                0014 C     | subroutine bling_light
                0015 C     | o calculate effective light for phytoplankton growth
                0016 C     |   There are multiple types of light.
                0017 C     | - irr_inst is the instantaneous irradiance field.
e0f9a7ba0b Matt*0018 C     | - irr_mix is the same, but with the irr_inst averaged throughout
                0019 C     |   the mixed layer. This quantity is intended to represent the
                0020 C     |   light to which phytoplankton subject to turbulent transport in
c0d1c06c15 Matt*0021 C     |   the mixed-layer would be exposed.
e0f9a7ba0b Matt*0022 C     | - irr_mem is a temporally smoothed field carried between
c0d1c06c15 Matt*0023 C     |   timesteps, to represent photoadaptation.
e0f9a7ba0b Matt*0024 C     | - irr_eff is the effective irradiance for photosynthesis,
c0d1c06c15 Matt*0025 C     |   given either by irr_inst or irr_mix, depending on model
                0026 C     |   options and location.
e0f9a7ba0b Matt*0027 C     | o instantaneous light is calculated either from
                0028 C     | - date and latitude, then exponentially attenuated down the
                0029 C     |   water column, or
                0030 C     | - short-wave radiation read from external forcing file,
                0031 C     |   attenuated down the water column, or
                0032 C     | - short-wave radiation distributed through the water column
                0033 C     |   according to SWFRAC routine
c0d1c06c15 Matt*0034 C     =================================================================
                0035
e0f9a7ba0b Matt*0036       IMPLICIT NONE
                0037
c0d1c06c15 Matt*0038 C     === Global variables ===
                0039 C     irr_inst      :: Instantaneous irradiance
                0040 C     irr_mem       :: Phyto irradiance memory
                0041
                0042 #include "SIZE.h"
e0f9a7ba0b Matt*0043 #include "DYNVARS.h"
c0d1c06c15 Matt*0044 #include "EEPARAMS.h"
                0045 #include "PARAMS.h"
                0046 #include "GRID.h"
                0047 #include "BLING_VARS.h"
e0f9a7ba0b Matt*0048 #ifdef USE_QSW
                0049 #include "FFIELDS.h"
                0050 #endif
a284455135 Matt*0051 #ifdef ALLOW_AUTODIFF_TAMC
c0d1c06c15 Matt*0052 # include "tamc.h"
                0053 #endif
                0054
                0055 C     === Routine arguments ===
                0056 C     bi,bj         :: tile indices
                0057 C     iMin,iMax     :: computation domain: 1rst index range
                0058 C     jMin,jMax     :: computation domain: 2nd  index range
                0059 C     myTime        :: current time
                0060 C     myIter        :: current timestep
                0061 C     myThid        :: thread Id. number
                0062       INTEGER bi, bj, imin, imax, jmin, jmax
                0063       INTEGER myThid
                0064       INTEGER myIter
                0065       _RL     myTime
                0066 C     === Input ===
                0067       _RL mld       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068 C     === Output ===
e0f9a7ba0b Matt*0069 C      irr_inst     :: instantaneous light
c0d1c06c15 Matt*0070 C      irr_eff      :: effective light for photosynthesis
                0071       _RL irr_inst  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0072       _RL irr_eff   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0073
                0074 C     === Local variables ===
                0075       _RL atten
                0076       _RL irr_surf  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e0f9a7ba0b Matt*0077 #ifdef ML_MEAN_LIGHT
c0d1c06c15 Matt*0078       _RL irr_mix   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL SumMLIrr
                0080       _RL tmp_ML
                0081 #endif
                0082 #ifndef USE_QSW
9f5240b52a Jean*0083       _RL solar, albedo
                0084       _RL dayfrac, yday, delta
                0085       _RL lat, sun1, dayhrs
                0086       _RL cosz, frac, fluxi
c0d1c06c15 Matt*0087       _RL sfac      (1-OLy:sNy+OLy)
                0088 #endif
00fa2d4ddd mmaz*0089 #ifdef PHYTO_SELF_SHADING
                0090       _RL k0_rd, chi_rd, e_rd
                0091       _RL k0_bg, chi_bg, e_bg
a284455135 Matt*0092       _RL kChl_rd   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0093       _RL kChl_bg   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
00fa2d4ddd mmaz*0094       _RL atten_rd
                0095       _RL atten_bg
a284455135 Matt*0096       _RL irr_rd    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0097       _RL irr_bg    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
4e4ad91a39 Jean*0098 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0099 C     tkey :: tape key (tile dependent)
                0100       INTEGER tkey
00fa2d4ddd mmaz*0101 #endif
4e4ad91a39 Jean*0102 #endif /* PHYTO_SELF_SHADING */
e0f9a7ba0b Matt*0103       INTEGER i,j,k
a284455135 Matt*0104       LOGICAL QSW_underice
c0d1c06c15 Matt*0105 CEOP
                0106
a284455135 Matt*0107 #ifdef PHYTO_SELF_SHADING
                0108 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0109       tkey = bi + (bj - 1)*nSx + (ikey_dynamics - 1)*nSx*nSy
a284455135 Matt*0110 # endif /* ALLOW_AUTODIFF_TAMC */
                0111 #endif /* PHYTO_SELF_SHADING */
                0112
c0d1c06c15 Matt*0113        DO k=1,Nr
                0114         DO j=jmin,jmax
                0115           DO i=imin,imax
a284455135 Matt*0116               irr_eff(i,j,k) = 0. _d 0
                0117 #ifdef PHYTO_SELF_SHADING
                0118               irr_rd(i,j,k)        = 0. _d 0
                0119               irr_bg(i,j,k)        = 0. _d 0
                0120 #endif
c0d1c06c15 Matt*0121           ENDDO
                0122         ENDDO
                0123        ENDDO
                0124
                0125 c ---------------------------------------------------------------------
                0126 c  Surface insolation
                0127
e0f9a7ba0b Matt*0128 #ifndef USE_QSW
c0d1c06c15 Matt*0129 c  From pkg/dic/dic_insol
                0130 c  find light as function of date and latitude
                0131 c  based on paltridge and parson
                0132
                0133       solar  = 1360. _d 0   !solar constant
                0134       albedo = 0.6 _d 0     !planetary albedo
                0135
                0136 C     Case where a 2-d output array is needed: for now, stop here.
                0137       IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
                0138        STOP 'ABNORMAL END: S/R INSOL: 2-D output not implemented'
                0139       ENDIF
                0140
                0141 C find day (****NOTE for year starting in winter*****)
                0142         dayfrac=mod(myTime,360. _d 0*86400. _d 0)
                0143      &                    /(360. _d 0*86400. _d 0)  !fraction of year
                0144         yday = 2. _d 0*PI*dayfrac                    !convert to radians
                0145         delta = (0.006918 _d 0
                0146      &         -(0.399912 _d 0*cos(yday))            !cosine zenith angle
                0147      &         +(0.070257 _d 0*sin(yday))            !(paltridge+platt)
                0148      &         -(0.006758 _d 0*cos(2. _d 0*yday))
                0149      &         +(0.000907 _d 0*sin(2. _d 0*yday))
                0150      &         -(0.002697 _d 0*cos(3. _d 0*yday))
                0151      &         +(0.001480 _d 0*sin(3. _d 0*yday)) )
                0152        DO j=1-OLy,sNy+OLy
                0153 C latitude in radians
                0154           lat=YC(1,j,1,bj)*deg2rad
                0155 C     latitute in radians, backed out from coriolis parameter
                0156 C     (makes latitude independent of grid)
                0157           IF ( usingCartesianGrid .OR. usingCylindricalGrid )
                0158      &         lat = asin( fCori(1,j,1,bj)/(2. _d 0*omega) )
                0159           sun1 = -sin(delta)/cos(delta) * sin(lat)/cos(lat)
                0160           IF (sun1.LE.-0.999 _d 0) sun1=-0.999 _d 0
                0161           IF (sun1.GE. 0.999 _d 0) sun1= 0.999 _d 0
                0162           dayhrs = abs(acos(sun1))
                0163           cosz = ( sin(delta)*sin(lat)+              !average zenith angle
                0164      &            (cos(delta)*cos(lat)*sin(dayhrs)/dayhrs) )
                0165           IF (cosz.LE.5. _d -3) cosz= 5. _d -3
                0166           frac = dayhrs/PI                           !fraction of daylight in day
                0167 C daily average photosynthetically active solar radiation just below surface
e0f9a7ba0b Matt*0168           fluxi = solar*(1. _d 0-albedo)*cosz*frac*parfrac
c0d1c06c15 Matt*0169
                0170 C convert to sfac
                0171           sfac(j) = MAX(1. _d -5,fluxi)
                0172        ENDDO !j
                0173
9f5240b52a Jean*0174 #endif /* ndef USE_QSW */
c0d1c06c15 Matt*0175
                0176 c ---------------------------------------------------------------------
                0177 c  instantaneous light, mixed layer averaged light
                0178
                0179       DO j=jmin,jmax
                0180        DO i=imin,imax
e0f9a7ba0b Matt*0181
c0d1c06c15 Matt*0182 c  Photosynthetically-available radiations (PAR)
e0f9a7ba0b Matt*0183 #ifdef USE_QSW
c0d1c06c15 Matt*0184         irr_surf(i,j) = max(epsln,
                0185      &                 -parfrac*Qsw(i,j,bi,bj)*maskC(i,j,1,bi,bj))
                0186 #else
                0187         irr_surf(i,j) = sfac(j)
                0188 #endif
e0f9a7ba0b Matt*0189
                0190 c  Remove light under ice
                0191 c  If using Qsw and seaice/thsice, then ice fraction is already
                0192 c  taken into account
                0193         QSW_underice = .FALSE.
                0194 #ifdef USE_QSW
                0195         IF ( useSEAICE ) QSW_underice = .TRUE.
                0196         IF ( useThSIce ) QSW_underice = .TRUE.
                0197 #endif
                0198         IF ( .NOT. QSW_underice ) THEN
                0199          irr_surf(i,j) = irr_surf(i,j)*(1. _d 0 - FIce(i,j,bi,bj))
                0200         ENDIF
c0d1c06c15 Matt*0201
                0202 #ifdef ML_MEAN_LIGHT
                0203         SumMLIrr   = 0. _d 0
                0204         tmp_ML     = 0. _d 0
                0205 #endif
                0206
                0207         DO k=1,Nr
                0208
a284455135 Matt*0209 #ifdef PHYTO_SELF_SHADING
                0210 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0211 CADJ STORE irr_bg = comlev1_bibj, key=tkey, kind=isbyte
                0212 CADJ STORE irr_rd = comlev1_bibj, key=tkey, kind=isbyte
a284455135 Matt*0213 # endif /* ALLOW_AUTODIFF_TAMC */
                0214 #endif /* PHYTO_SELF_SHADING */
e0f9a7ba0b Matt*0215
a284455135 Matt*0216          IF (hFacC(i,j,k,bi,bj).gt.0) THEN
00fa2d4ddd mmaz*0217 #ifdef PHYTO_SELF_SHADING
                0218
                0219 c  Use bio-optical model of Manizza et al. (2005) to account for
6ffd1aa797 Jean*0220 c  effect of self-shading on ligt available for phytoplankton
00fa2d4ddd mmaz*0221 c  growth. As written this DOES NOT feedback onto the absorption
6ffd1aa797 Jean*0222 c  of shortwave radiation calculated in the physical model, which
00fa2d4ddd mmaz*0223 c  is instead calculated in the subroutine swfrac
                0224
                0225 c  Specify co-efficients for bio-optical model {kChl = k0 +chi[chl]^e}
                0226 c  in red and blue-green fractions (Morel 1988; Foujols et al. 2000)
                0227         k0_rd  = 0.225 _d 0
                0228         k0_bg  = 0.0232 _d 0
                0229         chi_rd = 0.037 _d 0
                0230         chi_bg = 0.074 _d 0
                0231         e_rd   = 0.629 _d 0
                0232         e_bg   = 0.674 _d 0
                0233
                0234          IF (k.eq.1) THEN
                0235 c  Attenuation coefficient adjusted to chlorophyll in top layer
                0236           kChl_rd(i,j,1) = k0_rd + chi_rd*(chl(i,j,1,bi,bj)**e_rd)
                0237           kChl_bg(i,j,1) = k0_bg + chi_bg*(chl(i,j,1,bi,bj)**e_bg)
                0238 c  Light attenuation in middle of top layer
                0239           atten_rd = kChl_rd(i,j,1)*drF(1)/2. _d 0*hFacC(i,j,1,bi,bj)
                0240           atten_bg = kChl_bg(i,j,1)*drF(1)/2. _d 0*hFacC(i,j,1,bi,bj)
                0241 c  Irradiance in middle of top layer
a284455135 Matt*0242           irr_rd(i,j,1) = irr_surf(i,j) * exp(-atten_rd) * 0.5 _d 0
                0243           irr_bg(i,j,1) = irr_surf(i,j) * exp(-atten_bg) * 0.5 _d 0
                0244           irr_inst(i,j,1) = irr_rd(i,j,1) + irr_bg(i,j,1)
00fa2d4ddd mmaz*0245          ELSE
                0246
                0247 c  Attenuation coefficient adjusted to chlorophyll in kth layer
                0248           kChl_rd(i,j,k) = k0_rd + chi_rd*(chl(i,j,k,bi,bj)**e_rd)
                0249           kChl_bg(i,j,k) = k0_bg + chi_bg*(chl(i,j,k,bi,bj)**e_bg)
                0250 c  Light attenuation from one more layer
                0251           atten_rd = kChl_rd(i,j,k)*drF(k)/2. _d 0*hFacC(i,j,k,bi,bj)
                0252      &       + kChl_rd(i,j,k-1)*drF(k-1)/2. _d 0*hFacC(i,j,k-1,bi,bj)
                0253           atten_bg = kChl_bg(i,j,k)*drF(k)/2. _d 0*hFacC(i,j,k,bi,bj)
                0254      &       + kChl_bg(i,j,k-1)*drF(k-1)/2. _d 0*hFacC(i,j,k-1,bi,bj)
                0255 c  Irradiance in middle of layer k
                0256           irr_rd(i,j,k) = irr_rd(i,j,k-1)*exp(-atten_rd)
6ffd1aa797 Jean*0257           irr_bg(i,j,k) = irr_bg(i,j,k-1)*exp(-atten_bg)
a284455135 Matt*0258           irr_inst(i,j,k) = irr_rd(i,j,k) + irr_bg(i,j,k)
00fa2d4ddd mmaz*0259          ENDIF
                0260
                0261 #else
e0f9a7ba0b Matt*0262 C SW radiation attenuated exponentially
c0d1c06c15 Matt*0263          IF (k.eq.1) THEN
                0264 c  Light attenuation in middle of top layer
                0265           atten = k0*drF(1)/2. _d 0*hFacC(i,j,1,bi,bj)
                0266           irr_inst(i,j,1) = irr_surf(i,j)*exp(-atten)
                0267          ELSE
                0268 c  Attenuation from one more layer
                0269           atten = k0*drF(k)/2. _d 0*hFacC(i,j,k,bi,bj)
                0270      &           + k0*drF(k-1)/2. _d 0*hFacC(i,j,k-1,bi,bj)
                0271           irr_inst(i,j,k) =
                0272      &           irr_inst(i,j,k-1)*exp(-atten)
                0273          ENDIF
                0274
00fa2d4ddd mmaz*0275 #endif /* if BLING_USE_SHADING */
e0f9a7ba0b Matt*0276
c0d1c06c15 Matt*0277 #ifdef ML_MEAN_LIGHT
                0278 c  Mean irradiance in the mixed layer
                0279          IF ((-rf(k+1) .le. mld(i,j)).and.
00fa2d4ddd mmaz*0280      &               (-rf(k+1).lt.MLmix_max)) THEN
c0d1c06c15 Matt*0281           SumMLIrr = SumMLIrr+drF(k)*irr_inst(i,j,k)
                0282           tmp_ML = tmp_ML + drF(k)
                0283           irr_mix(i,j) = SumMLIrr/tmp_ML
                0284          ENDIF
                0285 #endif
                0286
                0287          ENDIF
                0288
                0289         ENDDO
                0290        ENDDO
                0291       ENDDO
                0292
                0293       DO k=1,Nr
                0294        DO j=jmin,jmax
e0f9a7ba0b Matt*0295         DO i=imin,imax
c0d1c06c15 Matt*0296
                0297          IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
                0298           irr_eff(i,j,k) = irr_inst(i,j,k)
e0f9a7ba0b Matt*0299
c0d1c06c15 Matt*0300 #ifdef ML_MEAN_LIGHT
e0f9a7ba0b Matt*0301 c  Inside mixed layer, effective light is set to mean mixed layer light
c0d1c06c15 Matt*0302          IF ((-rf(k+1) .le. mld(i,j)).and.
00fa2d4ddd mmaz*0303      &               (-rf(k+1).lt.MLmix_max)) THEN
c0d1c06c15 Matt*0304            irr_eff(i,j,k) = irr_mix(i,j)
                0305           ENDIF
e0f9a7ba0b Matt*0306 #endif
c0d1c06c15 Matt*0307
                0308          ENDIF
                0309
                0310         ENDDO
                0311        ENDDO
                0312       ENDDO
e0f9a7ba0b Matt*0313
c0d1c06c15 Matt*0314       RETURN
                0315       END