/[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.3 - (show annotations) (download)
Tue Apr 16 20:20:27 2013 UTC (12 years, 3 months ago) by jahn
Branch: MAIN
Changes since 1.2: +5 -5 lines
fix types in format strings

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

  ViewVC Help
Powered by ViewVC 1.1.22