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

Annotation of /MITgcm_contrib/darwin2/pkg/darwin/darwin_readparams.F

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


Revision 1.7 - (hide annotations) (download)
Tue Apr 16 20:21:21 2013 UTC (12 years, 3 months ago) by jahn
Branch: MAIN
Changes since 1.6: +2 -2 lines
compute missing waveband information from provided values

1 jahn 1.7 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/darwin/darwin_readparams.F,v 1.6 2013/04/16 20:20:16 jahn Exp $
2 jahn 1.2 C $Name: $
3 jahn 1.1
4     #include "DARWIN_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: DARWIN_READPARMS
8    
9     C !INTERFACE: ==========================================================
10     SUBROUTINE DARWIN_READPARMS( myThid )
11    
12     C !DESCRIPTION:
13     C Initialize DARWIN parameters, read in data.darwin
14    
15     C !USES: ===============================================================
16     IMPLICIT NONE
17     #include "SIZE.h"
18     #include "EEPARAMS.h"
19     #include "PARAMS.h"
20     #include "GCHEM.h"
21     #include "DARWIN_SIZE.h"
22     #include "DARWIN_IO.h"
23     #include "DARWIN_PARAMS.h"
24     #ifdef WAVEBANDS
25     #include "SPECTRAL_SIZE.h"
26     #include "SPECTRAL.h"
27     #include "WAVEBANDS_PARAMS.h"
28     #endif
29    
30     C !INPUT PARAMETERS: ===================================================
31     C myThid :: thread number
32     INTEGER myThid
33    
34     C !OUTPUT PARAMETERS: ==================================================
35     C none
36    
37     #ifdef ALLOW_DARWIN
38    
39     C !LOCAL VARIABLES: ====================================================
40     C iUnit :: unit number for I/O
41     C msgBuf :: message buffer
42     INTEGER iUnit,errIO
43     CHARACTER*(MAX_LEN_MBUF) msgBuf
44     #if defined(WAVEBANDS) || defined(OASIM)
45     INTEGER ilam,i
46     _RL planck, c, hc, oavo, hcoavo, rlamm
47     #endif
48     CEOP
49    
50     NAMELIST /DARWIN_FORCING/
51     & darwin_iceFile,
52     & darwin_ironFile,
53     & darwin_PARFile,
54     & darwin_nutWVelFile,
55     & darwin_PO4_relaxFile, darwin_NO3_relaxFile,
56     & darwin_FeT_relaxFile, darwin_Si_relaxFile,
57     & darwin_PO4_fluxFile, darwin_NO3_fluxFile,
58     & darwin_FeT_fluxFile, darwin_Si_fluxFile,
59     & darwin_waterabsorbFile, darwin_phytoabsorbFile,
60     & darwin_surfacespecFile, darwin_acdomFile,
61     & darwin_particleabsorbFile,
62     #ifdef OASIM
63     & darwin_oasim_edFile, darwin_oasim_esFile,
64     #endif
65     & darwin_relaxscale,
66     & darwin_ForcingPeriod, darwin_ForcingCycle,
67 jahn 1.2 & darwin_PARunits, darwin_W_to_uEins,
68 jahn 1.3 & darwin_ironUnits,
69     & darwin_dustFrac,
70 jahn 1.1 #ifdef ALLOW_PAR_DAY
71     & darwin_PARavPeriod,
72     #endif
73     & darwin_seed
74    
75     #if defined(WAVEBANDS) || defined(OASIM)
76     NAMELIST /DARWIN_SPECTRAL_PARM/
77     & darwin_waves
78     #ifdef WAVEBANDS
79 jahn 1.4 & ,darwin_wavebands
80 jahn 1.1 #ifdef DAR_CALC_ACDOM
81     & ,darwin_Sdom
82     #endif
83     #ifdef DAR_DIAG_ACDOM
84     & ,darwin_diag_acdom_ilam
85     #endif
86     #ifdef DAR_RADTRANS
87     & ,darwin_PAR_ilamLo
88     & ,darwin_PAR_ilamHi
89     & ,darwin_radmodThresh
90     & ,darwin_Dmax
91     & ,darwin_rmus
92     & ,darwin_rmuu
93     & ,darwin_bbw
94     & ,darwin_bbphy
95     & ,darwin_bbmin
96     & ,darwin_radtrans_kmax
97     & ,darwin_radtrans_niter
98     & ,darwin_part_size_P
99     #endif
100     #endif /* WAVEBANDS */
101     #endif /* WAVEBANDS || OASIM */
102    
103     #ifdef DAR_DIAG_CHL
104     NAMELIST /DARWIN_CHL/
105     & Geider_Bigalphachl, Geider_smallalphachl,
106     & Geider_Bigchl2cmax, Geider_smallchl2cmax,
107     & Geider_Bigchl2cmin,
108     & Doney_Bmin, Doney_Bmax, Doney_PARstar,
109     & Cloern_A, Cloern_B, Cloern_C, Cloern_chl2cmin
110     #endif
111    
112     #ifdef ALLOW_CARBON
113     NAMELIST /DIC_FORCING/
114     & DIC_windFile, DIC_atmospFile,
115     & dic_pCO2, dic_int1, dic_int2, dic_int3, dic_int4
116     #endif
117    
118     C Set defaults values for parameters in DARWIN_IO.h
119     darwin_iceFile=' '
120     darwin_ironFile=' '
121     darwin_PARFile=' '
122     darwin_nutWVelFile=' '
123     darwin_PO4_relaxFile=' '
124     darwin_NO3_relaxFile=' '
125     darwin_FeT_relaxFile=' '
126     darwin_Si_relaxFile=' '
127     darwin_PO4_fluxFile=' '
128     darwin_NO3_fluxFile=' '
129     darwin_FeT_fluxFile=' '
130     darwin_Si_fluxFile=' '
131     darwin_waterabsorbFile=' '
132     darwin_phytoabsorbFile=' '
133     darwin_particleabsorbFile=' '
134     darwin_surfacespecFile=' '
135     darwin_acdomFile=' '
136     darwin_oasim_edFile=' '
137     darwin_oasim_esFile=' '
138 jahn 1.2 darwin_PARunits='Ein/m2/d '
139     darwin_W_to_uEins=1. _d 0/0.2174 _d 0
140 jahn 1.3 darwin_ironUnits='mol Fe/m2/s '
141     darwin_dustFrac=0.035 _d 0
142 jahn 1.1 darwin_relaxscale=0. _d 0
143     darwin_seed=0
144     c default periodic forcing to same as for GCHEM
145     darwin_forcingPeriod=gchem_ForcingPeriod
146     darwin_forcingCycle=gchem_ForcingCycle
147    
148     #ifdef ALLOW_CARBON
149     DIC_windFile = ' '
150     DIC_atmospFile= ' '
151     dic_int1 = 0
152     dic_int2 = 0
153     dic_int3 = 0
154     dic_int4 = 0
155     dic_pCO2 = 0. _d 0
156     #endif
157    
158     #ifdef ALLOW_PAR_DAY
159     darwin_PARavPeriod=86400. _d 0
160     #endif
161    
162     #if defined(WAVEBANDS) || defined(OASIM)
163     DO ilam=1,tlam
164     darwin_waves(ilam) = 0
165     ENDDO
166     #endif
167    
168     #ifdef WAVEBANDS
169 jahn 1.4 DO ilam=1,tlam+1
170 jahn 1.7 darwin_wavebands(ilam) = -1
171 jahn 1.4 ENDDO
172 jahn 1.1 #ifdef DAR_CALC_ACDOM
173     darwin_Sdom = 0.014 _d 0
174     #endif
175     #ifdef DAR_DIAG_ACDOM
176     c value >= 100 will be converted to index in wavebands_init_fixed
177     darwin_diag_acdom_ilam = 450
178     #endif
179     #ifdef DAR_RADTRANS
180     darwin_PAR_ilamLo = 1
181     darwin_PAR_ilamHi = tlam
182     darwin_radmodThresh = 1 _d -4
183     darwin_Dmax = 500 _d 0
184     darwin_rmus = 1.0/0.83 _d 0
185     darwin_rmuu = 1.0/0.4 _d 0
186     darwin_bbmin = 0.0002 _d 0
187     darwin_bbw = 0.5 _d 0
188     do i=1,tnabp
189     darwin_bbphy(i) = 0 _d 0
190     enddo
191     darwin_radtrans_kmax = Nr
192 jahn 1.5 darwin_radtrans_niter = -2
193 jahn 1.1 darwin_part_size_P = 1 _d -15 ! mmol P per particle
194     #endif
195     #endif /* WAVEBANDS */
196    
197     C Open and read the data.darwin file
198     _BEGIN_MASTER(myThid)
199     WRITE(msgBuf,'(A)') ' DARWIN_READPARMS: opening data.darwin'
200     CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,
201     & SQUEEZE_RIGHT , 1)
202     CALL OPEN_COPY_DATA_FILE(
203     I 'data.darwin', 'DARWIN_READPARAMS',
204     O iUnit,
205     I myThid )
206     READ(UNIT=iUnit,NML=DARWIN_FORCING)
207     #ifdef ALLOW_PAR_DAY
208     darwin_PARnav = NINT(darwin_PARavPeriod*nsubtime/dTtracerLev(1))
209     #endif
210 jahn 1.2 C factor for conversion to uEin/m2/s
211     IF ( darwin_PARunits(1:16) .EQ. 'uEin/m2/s ' ) THEN
212     darwin_PARFileConv = 1. _d 0
213     ELSEIF ( darwin_PARunits(1:16) .EQ. 'Ein/m2/d ' ) THEN
214     darwin_PARFileConv = 1. _d 6/86400. _d 0
215     ELSEIF ( darwin_PARunits(1:16) .EQ. 'W/m2 ' ) THEN
216     darwin_PARFileConv = darwin_W_to_uEins
217     ELSE
218     WRITE(msgBuf,'(2A)') 'S/R DARWIN_READPARMS:',
219     & 'darwin_PARunits must be one of Ein/m2/d, uEin/m2/s, W/m2'
220     CALL PRINT_ERROR( msgBuf , 1)
221     STOP 'unknown darwin_PARunits'
222     ENDIF
223 jahn 1.1
224 jahn 1.3 C factor for conversion to mmol Fe/m2/s (bioavailable fraction is done later)
225     IF ( darwin_ironUnits(1:16) .EQ. 'mol Fe/m2/s ' ) THEN
226     darwin_ironFileConv = 1000. _d 0
227     ELSEIF ( darwin_ironUnits(1:16) .EQ. 'kg dust/m2/s ' ) THEN
228     darwin_ironFileConv = 1000.*darwin_dustFrac*1000./55.845 _d 0
229     ELSE
230     WRITE(msgBuf,'(2A)') 'S/R DARWIN_READPARMS:',
231     & 'darwin_ironUnits must be one of mol Fe/m2/s, kg dust/m2/s'
232     CALL PRINT_ERROR( msgBuf , 1)
233     STOP 'unknown darwin_Feunits'
234     ENDIF
235    
236 jahn 1.1 #ifdef DAR_DIAG_CHL
237     C default values
238     C Geider: chl:c = max(chl2cmin, chl2cmax/(1+(chl2cmax*alphachl*PARday)/(2*Pcm)))
239     C Pcm = mu*limit*phytoTempFunction
240     Geider_smallalphachl = 2. _d -6 ! mmol C (mg Chl)-1 m2 (uEin)-1
241     Geider_Bigalphachl = 1. _d -6 ! mmol C (mg Chl)-1 m2 (uEin)-1
242     Geider_smallchl2cmax = 0.35 _d 0 ! mg Chl (mmol C)-1
243     Geider_Bigchl2cmax = 0.65 _d 0 ! mg Chl (mmol C)-1
244     Geider_smallchl2cmin = 0.003 _d 0 * 12. _d 0 ! mg Chl a/mmol C
245     Geider_Bigchl2cmin = 0.003 _d 0 * 12. _d 0 ! mg Chl a/mmol C
246    
247     C Doney: chl:c = (Bmax - (Bmax-Bmin)*MIN(1,PARday/PARstar))*limit
248     Doney_Bmax = 12. _d 0 / 37. _d 0 ! mg Chl a/mmol C
249     Doney_Bmin = 12. _d 0 / 90. _d 0 ! mg Chl a/mmol C
250     Doney_PARstar = 90. _d 0 / 0.2174 _d 0 ! uEin/m2/s
251    
252     C Cloern: chl:c = chl2cmin + A*exp(B*T)*exp(-C*PARday)*limit
253     Cloern_chl2cmin = 0.003 _d 0 * 12. _d 0 ! mg Chl a/mmol C
254     Cloern_A = 0.0154 _d 0 * 12. _d 0 ! mg Chl a/mmol C
255     Cloern_B = 0.050 _d 0 ! (degree C)^{-1}
256     Cloern_C = 0.059 _d 0 * 86400. _d 0 / 1. _d 6 ! m^2 s/uEin
257    
258     READ(UNIT=iUnit,NML=DARWIN_CHL)
259     #endif /* DAR_DIAG_CHL */
260    
261     #ifdef ALLOW_CARBON
262     READ(UNIT=iUnit,NML=DIC_FORCING)
263     #endif
264    
265     #ifdef DAR_RADTRANS
266     #ifndef DAR_NONSPECTRAL_BACKSCATTERING_RATIO
267     DO i=1,tnabp
268     IF ( darwin_bbphy(i) .NE. 0 _d 0 ) THEN
269     WRITE(msgBuf,'(2A)') 'S/R DARWIN_READPARMS:',
270     & 'darwin_bbphy is obsolete.'
271     CALL PRINT_ERROR( msgBuf , 1)
272     WRITE(msgBuf,'(2A)') 'S/R DARWIN_READPARMS:',
273     & 'Backscattering coefficients are now read from'
274     CALL PRINT_ERROR( msgBuf , 1)
275     WRITE(msgBuf,'(2A)') 'S/R DARWIN_READPARMS:',
276     & 'darwin_phytoabsorbFile.'
277     CALL PRINT_ERROR( msgBuf , 1)
278     ENDIF
279     ENDDO
280     #endif
281     #endif
282    
283     #if defined(WAVEBANDS) || defined(OASIM)
284     READ(UNIT=iUnit,NML=DARWIN_SPECTRAL_PARM,IOSTAT=errIO)
285     IF ( errIO .LT. 0 ) THEN
286     WRITE(msgBuf,'(A)')
287     & 'S/R DARWIN_READPARMS'
288     CALL PRINT_ERROR( msgBuf , 1)
289     WRITE(msgBuf,'(A)')
290     & 'Error reading darwin package'
291     CALL PRINT_ERROR( msgBuf , 1)
292     WRITE(msgBuf,'(A)')
293     & 'parameter file "data.darwin"'
294     CALL PRINT_ERROR( msgBuf , 1)
295     WRITE(msgBuf,'(A)')
296     & 'Problem in namelist DARWIN_SPECTRAL_PARM'
297     CALL PRINT_ERROR( msgBuf , 1)
298     STOP 'ABNORMAL END: S/R DARWIN_READPARMS'
299     ENDIF
300     #endif
301    
302     WRITE(msgBuf,'(A)')
303     & ' DARWIN_READPARMS: finished reading data.darwin'
304     CALL PRINT_MESSAGE(msgBuf, standardMessageUnit,
305     & SQUEEZE_RIGHT , 1)
306    
307     C Close the open data file
308     CLOSE(iUnit)
309    
310     #if defined(WAVEBANDS) || defined(OASIM)
311 jahn 1.4 IF ( darwin_waves(1).EQ.0.0 .AND. tlam.EQ.13 ) THEN
312     darwin_waves(1) = 400
313     darwin_waves(2) = 425
314     darwin_waves(3) = 450
315     darwin_waves(4) = 475
316     darwin_waves(5) = 500
317     darwin_waves(6) = 525
318     darwin_waves(7) = 550
319     darwin_waves(8) = 575
320     darwin_waves(9) = 600
321     darwin_waves(10) = 625
322     darwin_waves(11) = 650
323     darwin_waves(12) = 675
324     darwin_waves(13) = 700
325     IF ( darwin_wavebands(1).NE.0.0 ) THEN
326     WRITE(msgBuf,'(2A)') 'DARWIN_READPARMS: ',
327     & 'darwin_wavebands given without darwin_waves.'
328     CALL PRINT_ERROR( msgBuf, myThid )
329     WRITE(msgBuf,'(2A)') 'DARWIN_READPARMS: ',
330     & 'using standard waveband centers 400, 425, ...'
331     CALL PRINT_ERROR( msgBuf, myThid )
332     WRITE(msgBuf,'(2A)') 'DARWIN_READPARMS: ',
333     & 'provide darwin_waves if this is not desired.'
334     CALL PRINT_ERROR( msgBuf, myThid )
335     ENDIF
336     ENDIF
337 jahn 1.1 c Quanta conversion
338     planck = 6.6256 _d -34 !Plancks constant J sec
339     c = 2.998 _d 8 !speed of light m/sec
340     hc = 1.0/(planck*c)
341     oavo = 1.0/6.023 _d 23 ! 1/Avogadros number
342     hcoavo = hc*oavo
343     do ilam = 1,tlam
344     rlamm = darwin_waves(ilam)*1 _d -9 !lambda in m
345     WtouEins(ilam) = 1 _d 6*rlamm*hcoavo !Watts to uEin/s conversion
346 jahn 1.4 IF ( rlamm .EQ. 0.0 ) THEN
347     WRITE(msgBuf,'(2A)') 'DARWIN_READPARMS: ',
348     & 'please provide wavelengths in darwin_waves.'
349     CALL PRINT_ERROR( msgBuf, myThid )
350     STOP 'ABNORMAL END: S/R DARWIN_READPARMS'
351     ENDIF
352 jahn 1.1 enddo
353     #endif
354    
355     C-- Print a summary of parameter values:
356     iUnit = standardMessageUnit
357     WRITE(msgBuf,'(A)') '// ==================================='
358     CALL PRINT_MESSAGE( msgBuf, iUnit, SQUEEZE_RIGHT , myThid )
359     WRITE(msgBuf,'(A)') '// darwin parameters '
360     CALL PRINT_MESSAGE( msgBuf, iUnit, SQUEEZE_RIGHT , myThid )
361     WRITE(msgBuf,'(A)') '// ==================================='
362     CALL PRINT_MESSAGE( msgBuf, iUnit, SQUEEZE_RIGHT , myThid )
363     CALL WRITE_0D_I( darwin_seed, INDEX_NONE,
364     & 'darwin_seed =',
365     & ' /* seed for random number generator */')
366    
367     WRITE(msgBuf,'(A)') ' -----------------------------------'
368     CALL PRINT_MESSAGE( msgBuf, iUnit, SQUEEZE_RIGHT, myThid )
369    
370     #if defined(WAVEBANDS) || defined(OASIM)
371 jahn 1.6 CALL WRITE_1D_RL( darwin_waves, tlam, 0,
372 jahn 1.1 & 'darwin_waves =',
373     &' /* "central" wavelengths of wavebands */')
374     #endif
375     #ifdef WAVEBANDS
376     #ifdef DAR_CALC_ACDOM
377     CALL WRITE_0D_RL( darwin_Sdom, INDEX_NONE,
378     & 'darwin_Sdom =',
379     &' /* spectral slope for aCDOM */')
380     #endif
381     #ifdef DAR_DIAG_ACDOM
382     CALL WRITE_0D_I( darwin_diag_acdom_ilam, INDEX_NONE,
383     & 'darwin_diag_acdom_ilam =',
384     &' /* waveband to use for aCDOM diagnostic */')
385     #endif
386     #ifdef DAR_RADTRANS
387     CALL WRITE_0D_I( darwin_PAR_ilamLo, INDEX_NONE,
388     & 'darwin_PAR_ilamLo =',
389     &' /* waveband index of PAR lower bound */')
390     CALL WRITE_0D_I( darwin_PAR_ilamHi, INDEX_NONE,
391     & 'darwin_PAR_ilamHi =',
392     &' /* waveband index of PAR upper bound */')
393     CALL WRITE_0D_RL( darwin_radmodThresh, INDEX_NONE,
394     & 'darwin_radmodThresh =',
395     &' /* threshold for calling radmod (W/m2/waveband) */')
396     CALL WRITE_0D_RL( darwin_Dmax, INDEX_NONE,
397     & 'darwin_Dmax =',
398     &' /* depth at which Ed is assumed zero */')
399     CALL WRITE_0D_RL( darwin_rmus, INDEX_NONE,
400     & 'darwin_rmus =',
401     &' /* inverse average cosine of downward diffuse irradiance */')
402     CALL WRITE_0D_RL( darwin_rmuu, INDEX_NONE,
403     & 'darwin_rmuu =',
404     &' /* inverse average cosine of upward diffuse irradiance */')
405     CALL WRITE_0D_RL( darwin_bbw, INDEX_NONE,
406     & 'darwin_bbw =',
407     &' /* backscattering to forward scattering ratio for water */')
408     CALL WRITE_0D_RL( darwin_bbmin, INDEX_NONE,
409     & 'darwin_bbmin =',
410     &' /* minimum backscattering coefficient (1/m) */')
411     CALL WRITE_1D_RL( darwin_bbphy, tnabp, 0,
412     & 'darwin_bbphy =',
413     &' /* backscattering to forward scattering ratio for phyto */')
414     CALL WRITE_0D_I( darwin_radtrans_kmax, INDEX_NONE,
415     & 'darwin_radtrans_kmax =',
416     &' /* deepest level in which to compute irradiances */')
417     CALL WRITE_0D_I( darwin_radtrans_niter, INDEX_NONE,
418     & 'darwin_radtrans_niter =',
419 jahn 1.5 &' /* iterations/method for solveing 3-stream equations */')
420 jahn 1.1 #endif /* DAR_RADTRANS */
421     #endif /* WAVEBANDS */
422    
423     WRITE(msgBuf,'(A)') ' ==================================='
424     CALL PRINT_MESSAGE( msgBuf, iUnit, SQUEEZE_RIGHT, myThid )
425    
426     _END_MASTER(myThid)
427    
428     C Everyone else must wait for the parameters to be loaded
429     _BARRIER
430    
431     #endif /* ALLOW_DARWIN */
432    
433     RETURN
434     END

  ViewVC Help
Powered by ViewVC 1.1.22