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

Annotation 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 - (hide 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 jahn 1.3 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 jahn 1.2 C $Name: $
3 jahn 1.1
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 jahn 1.2 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 jahn 1.1 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 jahn 1.3 WRITE(msgBuf,'(2A,I3,A,I4,A,I4)') 'WAVEBANDS_INIT_FIXED: ',
104 jahn 1.1 & '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 jahn 1.3 WRITE(msgBuf,'(2A,I3,A,I4,A,I4)') 'WAVEBANDS_INIT_FIXED: ',
145 jahn 1.1 & '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 jahn 1.3 WRITE(msgBuf,'(2A,I3,A,I4,A,I4)') 'WAVEBANDS_INIT_FIXED: ',
182 jahn 1.1 & '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 jahn 1.3 WRITE(msgBuf,'(2A,I3,A,I4,A,I4)') 'WAVEBANDS_INIT_FIXED: ',
281 jahn 1.1 & '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