/[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.1 - (hide annotations) (download)
Wed Apr 13 18:56:25 2011 UTC (14 years, 3 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 jahn 1.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