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

Contents of /MITgcm_contrib/darwin2/pkg/darwin/wavebands_init_fixed.F

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


Revision 1.1 - (show annotations) (download)
Wed Apr 13 18:56:25 2011 UTC (14 years, 4 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt63c_20111011, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63p_20120707, ctrb_darwin2_ckpt63d_20111107, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_baseline, ctrb_darwin2_ckpt63n_20120604, ctrb_darwin2_ckpt63k_20120317, ctrb_darwin2_ckpt62z_20110622
darwin2 initial checkin

1 C $Header$
2 C $Name$
3
4 c ANNA wavebands_init_fixed.F reads-in and assigns input paramters for WAVEBANDS.
5
6 #include "DARWIN_OPTIONS.h"
7
8 CBOP
9 C !ROUTINE: WAVEBANDS_INIT_FIXED
10 C !INTERFACE:
11 subroutine wavebands_init_fixed(myThid)
12
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | SUBROUTINE WAVEBANDS_INIT_FIXED
16 C | o reads-in and assigns input paramters for WAVEBANDS.
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 implicit none
22 C == Global variables ===
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "DARWIN_SIZE.h"
27 #include "SPECTRAL_SIZE.h"
28 #include "SPECTRAL.h"
29 #include "WAVEBANDS_PARAMS.h"
30 #include "DARWIN_IO.h"
31
32 C !INPUT/OUTPUT PARAMETERS:
33 C == Routine arguments ==
34 C myThid :: my Thread Id number
35 integer myThid
36 CEOP
37
38 #ifdef WAVEBANDS
39
40 C !LOCAL VARIABLES:
41 C == Local variables ==
42 c local variables
43 CHARACTER*(MAX_LEN_MBUF) msgBuf
44 character*80 title
45 integer iUnit
46 integer swlambda,splambda,ssflambda
47 _RL sap,sap_ps,sbp,sbbp
48 _RL saw,sbw
49 _RL ssf
50 c _RL planck, c, hc, oavo, hcoavo, rlamm
51 #ifdef DAR_CALC_ACDOM
52 _RL rlam450,rlam
53 #else
54 _RL sacdom
55 #endif
56 c local indeces
57 integer nabp,i,ilam
58
59 do i = 1,tlam
60 pwaves(i) = darwin_waves(i)
61 enddo
62
63 C band widths used to convert OASIM data to irradiation per nm
64 C put boundaries half-way between central values
65 C but first and last boundary are at first and last "central" value
66 wb_width(1) = .5*(pwaves(2)-pwaves(1))
67 do i=2,tlam-1
68 wb_width(i) = .5*(pwaves(i+1)-pwaves(i-1))
69 enddo
70 wb_width(tlam) = .5*(pwaves(tlam)-pwaves(tlam-1))
71 wb_totalWidth = 0.0
72 do i=1,tlam
73 wb_totalWidth = wb_totalWidth + wb_width(i)
74 enddo
75 if (wb_totalWidth.LE.0) then
76 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: ',
77 & 'please provide wavelengths in darwin_waves.'
78 CALL PRINT_ERROR( msgBuf, myThid )
79 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
80 endif
81
82
83 c Water data files
84 if (darwin_waterabsorbFile .NE. ' ' ) THEN
85 CALL MDSFINDUNIT( iUnit, myThid )
86 open(iUnit,file=darwin_waterabsorbFile,
87 & status='old',form='formatted')
88 do i = 1,6 ! six lines of text for the header
89 read(iUnit,'(a50)')title ! trucates or pads (with spaces) to 50 characters length
90 enddo
91 do ilam = 1,tlam
92 read(iUnit,20)swlambda,saw,sbw
93 if (swlambda.NE.pwaves(ilam)) then
94 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: ',
95 & "wavelength for water spectrum doesn't match darwin_waves:"
96 CALL PRINT_ERROR( msgBuf, myThid )
97 WRITE(msgBuf,'(2A,I3,A,I4,A,F4.0)') 'WAVEBANDS_INIT_FIXED: ',
98 & 'ilam', ilam, ': ', swlambda, ' versus ', pwaves(ilam)
99 CALL PRINT_ERROR( msgBuf, myThid )
100 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
101 endif
102 aw(ilam) = saw
103 bw(ilam) = sbw
104 enddo
105 close(iUnit)
106 20 format(i5,f15.4,f10.4)
107 else
108 WRITE(msgBuf,'(A)')
109 & 'WAVEBANDS_INIT_FIXED: need to specify water absorption'
110 CALL PRINT_ERROR( msgBuf, myThid )
111 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
112 endif
113
114
115 c phyto data files
116 c ANNA phyto input data files must have a column for absorption by PS pigs
117 c ANNA easiest way to 'turn off' PS for growth is to put same values in both abs columns
118 if (darwin_phytoabsorbFile.NE. ' ' ) THEN
119 CALL MDSFINDUNIT( iUnit, myThid )
120 open(iUnit,file=darwin_phytoabsorbFile,
121 & status='old',form='formatted')
122 do i = 1,6 ! six lines of text for the header
123 read(iUnit,'(a50)')title
124 enddo
125 sbbp = 0. _d 0
126 do nabp = 1,tnabp
127 read(iUnit,'(a50)')title ! reads one line of text for the phytoplankton type header
128 do ilam = 1,tlam
129 #ifdef DAR_NONSPECTRAL_BACKSCATTERING_RATIO
130 read(iUnit,30)splambda,sap,sap_ps,sbp
131 #else
132 read(iUnit,'(i4,3f10.0,f20.0)')splambda,sap,sap_ps,sbp,sbbp
133 #endif
134 if (splambda.NE.pwaves(ilam)) then
135 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: ',
136 & "wavelength for phyto spectrum doesn't match darwin_waves:"
137 CALL PRINT_ERROR( msgBuf, myThid )
138 WRITE(msgBuf,'(2A,I3,A,I4,A,F4.0)') 'WAVEBANDS_INIT_FIXED: ',
139 & 'ilam', ilam, ': ', splambda, ' versus ', pwaves(ilam)
140 CALL PRINT_ERROR( msgBuf, myThid )
141 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
142 endif
143 ap(nabp,ilam) = sap
144 ap_ps(nabp,ilam) = sap_ps
145 bp(nabp,ilam) = sbp
146 bbp(nabp,ilam) = sbbp
147 enddo
148 enddo
149 close(iUnit)
150 30 format(i4,3f10.4)
151 else
152 WRITE(msgBuf,'(A)')
153 & 'WAVEBANDS_INIT_FIXED: need to specify phyto absorption'
154 CALL PRINT_ERROR( msgBuf, myThid )
155 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
156 endif
157
158
159 #ifndef OASIM
160 c QQ NEED IN HERE ifndef OASIM
161 c surface spectrum for initial use
162 if (darwin_surfacespecFile .NE. ' ' ) THEN
163 CALL MDSFINDUNIT( iUnit, myThid )
164 open(iUnit,file=darwin_surfacespecFile,
165 & status='old',form='formatted')
166 do i = 1,3 ! three lines of text for the header
167 read(iUnit,'(a50)')title
168 enddo
169 do ilam = 1,tlam
170 read(iUnit,40)ssflambda,ssf
171 if (ssflambda.NE.pwaves(ilam)) then
172 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: ',
173 & "wavelength for surface spectrum doesn't match darwin_waves:"
174 CALL PRINT_ERROR( msgBuf, myThid )
175 WRITE(msgBuf,'(2A,I3,A,I4,A,F4.0)') 'WAVEBANDS_INIT_FIXED: ',
176 & 'ilam', ilam, ': ', ssflambda, ' versus ', pwaves(ilam)
177 CALL PRINT_ERROR( msgBuf, myThid )
178 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
179 endif
180 sf(ilam) = ssf
181 enddo
182 close(iUnit)
183 40 format(i5,f15.6)
184 else
185 WRITE(msgBuf,'(A)')
186 & 'WAVEBANDS_INIT_FIXED: need surface spectrum'
187 CALL PRINT_ERROR( msgBuf, myThid )
188 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
189 endif
190 #endif /* not OASIM */
191
192
193 c absorption by cdom
194 #ifndef DAR_CALC_ACDOM
195 c if no file given then CDOM is zero
196 if (darwin_acdomFile.NE. ' ' ) THEN
197 CALL MDSFINDUNIT( iUnit, myThid )
198 open(iUnit,file=darwin_acdomFile,
199 & status='old',form='formatted')
200 do i = 1,6 ! six lines of text for the header
201 read(iUnit,'(a50)')title
202 enddo
203 do i = 1,tlam
204 read(iUnit,50)sacdom
205 acdom(i) = sacdom
206 enddo
207 close(iUnit)
208 50 format(f10.4)
209 else
210 WRITE(msgBuf,'(A)')
211 & 'WAVEBANDS_INIT_FIXED: no aCDOM'
212 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
213 & SQUEEZE_RIGHT, 1 )
214
215 do i = 1,tlam
216 acdom(i) = 0. _d 0
217 enddo
218 endif
219 #else /* DAR_CALC_ACDOM */
220 c for 3-D or for direct comparison to RADTRANS would need the same formulation for CDOM as in radtrans.
221 c CDOM absorption exponent
222 rlam450 = 450.0 _d 0
223 do ilam = 1,tlam
224 if (pwaves(ilam) .eq. 450) nl450 = ilam
225 rlam = float(pwaves(ilam))
226 excdom(ilam) = exp(-darwin_Sdom*(rlam-rlam450))
227 enddo
228
229 WRITE(msgBuf,'(A,1P1E20.12)')
230 & 'WAVEBANDS_INIT_FIXED: darwin_Sdom = ', darwin_Sdom
231 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
232 & SQUEEZE_RIGHT, 1 )
233 WRITE(msgBuf,'(A,I3)')
234 & 'WAVEBANDS_INIT_FIXED: nl450 = ', nl450
235 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
236 & SQUEEZE_RIGHT, 1 )
237 #endif /* DAR_CALC_ACDOM */
238
239 #ifdef DAR_DIAG_ACDOM
240 c find waveband index for acdom diagnostic
241 if (darwin_diag_acdom_ilam.GE.100) then
242 do ilam = 1,tlam
243 if (pwaves(ilam) .eq. darwin_diag_acdom_ilam) then
244 darwin_diag_acdom_ilam = ilam
245 goto 60
246 endif
247 enddo
248 60 continue
249 endif
250
251 WRITE(msgBuf,'(A,I3,A,I4)')
252 & 'WAVEBANDS_INIT_FIXED: aCDOM diag ilam = ',
253 & darwin_diag_acdom_ilam, ', lambda = ',
254 & pwaves(darwin_diag_acdom_ilam)
255 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
256 & SQUEEZE_RIGHT, 1 )
257 #endif
258
259 #ifdef DAR_RADTRANS
260 c absorption and scattering by particles
261 if (darwin_particleabsorbFile .NE. ' ' ) THEN
262 CALL MDSFINDUNIT( iUnit, myThid )
263 open(iUnit,file=darwin_particleabsorbFile,
264 & status='old',form='formatted')
265 do i = 1,6 ! six lines of text for the header
266 read(iUnit,'(a50)')title ! trucates or pads (with spaces) to 50 characters length
267 enddo
268 do ilam = 1,tlam
269 read(iUnit,'(I4,3F15.0)')splambda,sap,sbp,sbbp
270 if (splambda.NE.pwaves(ilam)) then
271 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: ',
272 & "wavelength for particle spectrum doesn't match darwin_waves:"
273 CALL PRINT_ERROR( msgBuf, myThid )
274 WRITE(msgBuf,'(2A,I3,A,I4,A,F4.0)') 'WAVEBANDS_INIT_FIXED: ',
275 & 'ilam', ilam, ': ', splambda, ' versus ', pwaves(ilam)
276 CALL PRINT_ERROR( msgBuf, myThid )
277 STOP 'ABNORMAL END: S/R WAVEBANDS_INIT_FIXED'
278 endif
279 apart(ilam) = sap
280 bpart(ilam) = sbp
281 bbpart(ilam) = sbbp
282 apart_P(ilam) = sap/darwin_part_size_P
283 bpart_P(ilam) = sbp/darwin_part_size_P
284 bbpart_P(ilam) = sbbp/darwin_part_size_P
285 enddo
286 close(iUnit)
287 else
288 do ilam = 1,tlam
289 apart(ilam) = 0. _d 0
290 bpart(ilam) = 0. _d 0
291 bbpart(ilam) = 0. _d 0
292 apart_P(ilam) = 0. _d 0
293 bpart_P(ilam) = 0. _d 0
294 bbpart_P(ilam) = 0. _d 0
295 enddo
296 endif
297
298 c print a summary
299 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: waveband widths:'
300 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
301 & SQUEEZE_RIGHT, 1 )
302 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
303 & ' lam width'
304 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
305 & SQUEEZE_RIGHT, 1 )
306 do ilam = 1,tlam
307 WRITE(msgBuf,'(A,I4,F15.6)') 'WAVEBANDS_INIT_FIXED: ',
308 & pwaves(ilam), wb_width(ilam)
309 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
310 & SQUEEZE_RIGHT, 1 )
311 enddo
312 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
313 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
314 & SQUEEZE_RIGHT, 1 )
315 c
316 #ifndef OASIM
317 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: surface spectrum:'
318 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
319 & SQUEEZE_RIGHT, 1 )
320 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
321 & ' lam sf'
322 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
323 & SQUEEZE_RIGHT, 1 )
324 do ilam = 1,tlam
325 WRITE(msgBuf,'(A,I4,F15.6)') 'WAVEBANDS_INIT_FIXED: ',
326 & pwaves(ilam), sf(ilam)
327 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
328 & SQUEEZE_RIGHT, 1 )
329 enddo
330 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
331 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
332 & SQUEEZE_RIGHT, 1 )
333 #endif
334 c
335 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: water spectra:'
336 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
337 & SQUEEZE_RIGHT, 1 )
338 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
339 & ' lam aw bw'
340 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
341 & SQUEEZE_RIGHT, 1 )
342 do ilam = 1,tlam
343 WRITE(msgBuf,'(A,I4,F15.4,F10.4)') 'WAVEBANDS_INIT_FIXED: ',
344 & pwaves(ilam), aw(ilam), bw(ilam)
345 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
346 & SQUEEZE_RIGHT, 1 )
347 enddo
348 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
349 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
350 & SQUEEZE_RIGHT, 1 )
351 c
352 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: phyto spectra:'
353 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
354 & SQUEEZE_RIGHT, 1 )
355 do nabp = 1,tnabp
356 WRITE(msgBuf,'(A,I4)') 'WAVEBANDS_INIT_FIXED: type ',nabp
357 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
358 & SQUEEZE_RIGHT, 1 )
359 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
360 & ' lam ap ap_ps bp bbp'
361 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
362 & SQUEEZE_RIGHT, 1 )
363 do ilam = 1,tlam
364 WRITE(msgBuf,'(A,I4,3F10.4,F20.9)') 'WAVEBANDS_INIT_FIXED: ',
365 & pwaves(ilam), ap(nabp,ilam), ap_ps(nabp,ilam),
366 & bp(nabp,ilam), bbp(nabp,ilam)
367 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
368 & SQUEEZE_RIGHT, 1 )
369 enddo
370 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
371 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
372 & SQUEEZE_RIGHT, 1 )
373 enddo
374 c
375 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: particulate spectra:'
376 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
377 & SQUEEZE_RIGHT, 1 )
378 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
379 & ' lam apart bpart bbpart'
380 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
381 & SQUEEZE_RIGHT, 1 )
382 do ilam = 1,tlam
383 WRITE(msgBuf,'(A,I4,1P3G15.6)')'WAVEBANDS_INIT_FIXED: ',
384 & pwaves(ilam), apart(ilam), bpart(ilam), bbpart(ilam)
385 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
386 & SQUEEZE_RIGHT, 1 )
387 enddo
388 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
389 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
390 & SQUEEZE_RIGHT, 1 )
391 c
392 WRITE(msgBuf,'(2A)') 'WAVEBANDS_INIT_FIXED: particulate spectra ',
393 & 'in phosphorus units:'
394 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
395 & SQUEEZE_RIGHT, 1 )
396 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
397 & ' lam apart_P bpart_P bbpart_P'
398 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
399 & SQUEEZE_RIGHT, 1 )
400 do ilam = 1,tlam
401 WRITE(msgBuf,'(A,I4,2F15.9,F15.12)') 'WAVEBANDS_INIT_FIXED: ',
402 & pwaves(ilam), apart_P(ilam), bpart_P(ilam), bbpart_P(ilam)
403 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
404 & SQUEEZE_RIGHT, 1 )
405 enddo
406 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
407 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
408 & SQUEEZE_RIGHT, 1 )
409 c
410 #ifndef DAR_CALC_ACDOM
411 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED: CDOM spectrum:'
412 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
413 & SQUEEZE_RIGHT, 1 )
414 WRITE(msgBuf,'(A,A)') 'WAVEBANDS_INIT_FIXED: ',
415 & ' lam aCDOM'
416 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
417 & SQUEEZE_RIGHT, 1 )
418 do ilam = 1,tlam
419 WRITE(msgBuf,'(A,I4,F10.4)') 'WAVEBANDS_INIT_FIXED: ',
420 & pwaves(ilam), acdom(ilam)
421 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
422 & SQUEEZE_RIGHT, 1 )
423 enddo
424 WRITE(msgBuf,'(A)') 'WAVEBANDS_INIT_FIXED:'
425 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
426 & SQUEEZE_RIGHT, 1 )
427 #endif
428
429 c constants
430 pid = DACOS(-1.0D0)
431 rad = 180.0D0/pid
432 #endif
433
434 #endif /* WAVEBANDS */
435
436 return
437 end
438

  ViewVC Help
Powered by ViewVC 1.1.22