/[MITgcm]/MITgcm_contrib/bling/pkg/bling_light.F
ViewVC logotype

Contents of /MITgcm_contrib/bling/pkg/bling_light.F

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


Revision 1.1 - (show annotations) (download)
Fri May 23 17:33:43 2014 UTC (11 years, 2 months ago) by mmazloff
Branch: MAIN
Adding package BLING

1 C $Header: $
2 C $Name: $
3
4 #include "BLING_OPTIONS.h"
5
6 CBOP
7 subroutine BLING_LIGHT(
8 U irr_eff,
9 I bi, bj, imin, imax, jmin, jmax,
10 I myIter, myTime, myThid )
11 C =================================================================
12 C | subroutine bling_light
13 C | o calculate effective light for phytoplankton growth
14 C | There are multiple types of light.
15 C | - irr_inst is the instantaneous irradiance field.
16 C | - irr_mix is the same, but with the irr_inst averaged throughout
17 C | the mixed layer. This quantity is intended to represent the
18 C | light to which phytoplankton subject to turbulent transport in
19 C | the mixed-layer would be exposed.
20 C | - irr_mem is a temporally smoothed field carried between
21 C | timesteps, to represent photoadaptation.
22 C | - irr_eff is the effective irradiance for photosynthesis,
23 C | given either by irr_inst or irr_mix, depending on model
24 C | options and location.
25 C =================================================================
26
27 implicit none
28
29 C === Global variables ===
30 C irr_inst :: Instantaneous irradiance
31 C irr_mem :: Phyto irradiance memory
32
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "FFIELDS.h"
37 #include "GRID.h"
38 #include "DYNVARS.h"
39 #include "BLING_VARS.h"
40 #include "PTRACERS_SIZE.h"
41 #include "PTRACERS_PARAMS.h"
42 #ifdef ALLOW_AUTODIFF_TAMC
43 # include "tamc.h"
44 #endif
45
46 C === Routine arguments ===
47 C bi,bj :: tile indices
48 C iMin,iMax :: computation domain: 1rst index range
49 C jMin,jMax :: computation domain: 2nd index range
50 C myTime :: current time
51 C myIter :: current timestep
52 C myThid :: thread Id. number
53 INTEGER bi, bj, imin, imax, jmin, jmax
54 INTEGER myThid
55 INTEGER myIter
56 _RL myTime
57 C === Output ===
58 C irr_eff :: effective light for photosynthesis
59 _RL irr_eff (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
60
61 C === Local variables ===
62 _RL solar, albedo
63 _RL dayfrac, yday, delta
64 _RL lat, sun1, dayhrs
65 _RL cosz, frac, fluxi
66 _RL atten
67 _RL irr_surf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 #ifdef ML_MEAN_LIGHT
69 _RL irr_mix (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RL SumMLIrr
71 _RL SumMLDepth
72 _RL dens_surf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL dens_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL delta_dens(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75 #endif
76 #ifndef READ_PAR
77 #ifndef USE_QSW
78 _RL sfac (1-OLy:sNy+OLy)
79 #endif
80 #endif
81 integer i,j,k
82 CEOP
83
84 #ifdef ML_MEAN_LIGHT
85 c ---------------------------------------------------------------------
86 c Mixed layer depth
87
88 c Surface density
89 CALL FIND_RHO_2D(
90 I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
91 I theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
92 O dens_surf,
93 I 1, bi, bj, myThid )
94
95 DO k=1,Nr
96 DO j=jmin,jmax
97 DO i=imin,imax
98 if (k.eq.1) then
99 delta_dens(i,j,1) = 0. _d 0
100 else
101 delta_dens(i,j,k) = 9999. _d 0
102 endif
103 ENDDO
104 ENDDO
105 ENDDO
106
107 DO k = 2,Nr
108
109 c Potential density
110 CALL FIND_RHO_2D(
111 I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
112 I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
113 O dens_z,
114 I k, bi, bj, myThid )
115
116 DO j=jmin,jmax
117 DO i=imin,imax
118 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
119 delta_dens(i,j,k) = dens_z(i,j)-dens_surf(i,j)
120 ENDIF
121 ENDDO
122 ENDDO
123 ENDDO
124 #endif
125
126 c ---------------------------------------------------------------------
127 c Surface insolation
128
129 #ifndef USE_EXFQSW
130 c From pkg/dic/dic_insol
131 c find light as function of date and latitude
132 c based on paltridge and parson
133
134 solar = 1360. _d 0 !solar constant
135 albedo = 0.6 _d 0 !planetary albedo
136
137 C Case where a 2-d output array is needed: for now, stop here.
138 IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
139 STOP 'ABNORMAL END: S/R INSOL: 2-D output not implemented'
140 ENDIF
141
142 C find day (****NOTE for year starting in winter*****)
143 dayfrac=mod(myTime,360. _d 0*86400. _d 0)
144 & /(360. _d 0*86400. _d 0) !fraction of year
145 yday = 2. _d 0*PI*dayfrac !convert to radians
146 delta = (0.006918 _d 0
147 & -(0.399912 _d 0*cos(yday)) !cosine zenith angle
148 & +(0.070257 _d 0*sin(yday)) !(paltridge+platt)
149 & -(0.006758 _d 0*cos(2. _d 0*yday))
150 & +(0.000907 _d 0*sin(2. _d 0*yday))
151 & -(0.002697 _d 0*cos(3. _d 0*yday))
152 & +(0.001480 _d 0*sin(3. _d 0*yday)) )
153 DO j=1-OLy,sNy+OLy
154 C latitude in radians
155 lat=YC(1,j,1,bj)*deg2rad
156 C latitute in radians, backed out from coriolis parameter
157 C (makes latitude independent of grid)
158 IF ( usingCartesianGrid .OR. usingCylindricalGrid )
159 & lat = asin( fCori(1,j,1,bj)/(2. _d 0*omega) )
160 sun1 = -sin(delta)/cos(delta) * sin(lat)/cos(lat)
161 IF (sun1.LE.-0.999 _d 0) sun1=-0.999 _d 0
162 IF (sun1.GE. 0.999 _d 0) sun1= 0.999 _d 0
163 dayhrs = abs(acos(sun1))
164 cosz = ( sin(delta)*sin(lat)+ !average zenith angle
165 & (cos(delta)*cos(lat)*sin(dayhrs)/dayhrs) )
166 IF (cosz.LE.5. _d -3) cosz= 5. _d -3
167 frac = dayhrs/PI !fraction of daylight in day
168 C daily average photosynthetically active solar radiation just below surface
169 fluxi = solar*(1. _d 0-albedo)*cosz*frac*parfrac
170
171 C convert to sfac
172 sfac(j) = MAX(1. _d -5,fluxi)
173 ENDDO !j
174
175 #endif
176
177 c ---------------------------------------------------------------------
178 c instantaneous light, mixed layer averaged light
179
180 C$TAF LOOP = parallel
181 DO j=jmin,jmax
182 C$TAF LOOP = parallel
183 DO i=imin,imax
184
185 c Photosynthetically-available radiations (PAR)
186 #ifdef USE_EXFQSW
187 irr_surf(i,j) = max(0. _d 0,
188 & -parfrac*Qsw(i,j,bi,bj)*maskC(i,j,1,bi,bj))
189 #else
190 irr_surf(i,j) = sfac(j)
191 #endif
192 IF ( .NOT. QSW_underice ) THEN
193 c if using Qsw but not seaice/thsice or coupled, then
194 c ice fraction needs to be taken into account
195 irr_surf(i,j) = irr_surf(i,j)*(1. _d 0 - FIce(i,j,bi,bj))
196 ENDIF
197
198 #ifdef ML_MEAN_LIGHT
199 SumMLIrr = 0. _d 0
200 SumMLDepth = 0. _d 0
201 #endif
202
203 c C$TAF init ml_stuff = static, Nr
204 DO k=1,Nr
205 c C$TAF STORE SumMLDepth = ml_stuff
206
207 IF (hFacC(i,j,k,bi,bj).gt.0) THEN
208
209 IF (k.eq.1) THEN
210 c Light attenuation in middle of top layer
211 atten = k0*drF(1)/2. _d 0*hFacC(i,j,1,bi,bj)
212 irr_inst(i,j,1,bi,bj) = irr_surf(i,j)*exp(-atten)
213 ELSE
214 c Attenuation from one more layer
215 atten = k0*drF(k)/2. _d 0*hFacC(i,j,k,bi,bj)
216 & + k0*drF(k-1)/2. _d 0*hFacC(i,j,k-1,bi,bj)
217 irr_inst(i,j,k,bi,bj) =
218 & irr_inst(i,j,k-1,bi,bj)*exp(-atten)
219 ENDIF
220
221 #ifdef ML_MEAN_LIGHT
222 c Mean irradiance in the mixed layer
223 IF (delta_dens(i,j,k) .LT. 0.03 _d 0) then
224 SumMLIrr = SumMLIrr+drF(k)*irr_inst(i,j,k,bi,bj)
225 SumMLDepth = SumMLDepth+drF(k)
226 irr_mix(i,j) = SumMLIrr/SumMLDepth
227 ENDIF
228 #endif
229
230 ENDIF
231
232 ENDDO
233 ENDDO
234 ENDDO
235
236 c ---------------------------------------------------------------------
237 C Phytoplankton photoadaptation to local light level
238 DO k=1,Nr
239 DO j=jmin,jmax
240 DO i=imin,imax
241
242 irr_eff(i,j,k) = irr_inst(i,j,k,bi,bj)
243 #ifdef ML_MEAN_LIGHT
244 c Inside mixed layer, effective light is set to mean mixed layer light
245 IF (delta_dens(i,j,k) .LT. 0.03 _d 0) THEN
246 irr_eff(i,j,k) = irr_mix(i,j)
247 ENDIF
248 #endif
249
250 irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
251 & (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
252 & min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
253
254 ENDDO
255 ENDDO
256 ENDDO
257
258 RETURN
259 END

  ViewVC Help
Powered by ViewVC 1.1.22