Appendix to the report to SNF entitled „The Significance of coccolithophore blooms fort he oceanic carbon cycle: A numerical modelling experiment simulating Emiliania huxleyi blooms in the North Atlantic“

By Michael Knappertsbusch, developed from April 1, 1993 through September 30, 1994.
Date of report: November 2, 1994
Grant No. 8220-033288





Program source code to COCDIA 10.04, input data files to COCDIA, and source codes to auxiliary programs PP400, REDUC2, and TEST



Contents

Source code for program COCDIA, version 10.04					Page  2
Include file COCDAT.inc (valid for COCDIA and PP programs)			Page 27
Alternative form of function GMAXE(TEMP)					Page 29
Input files to the COCDIA program						Page 30
        Input file INITCON.DAT with initial conditions			Page 30
	Input file MIXED_LAYER_DEPTHS						Page 30
	Input file SEA_SURFACE_ALBEDO						Page 31
	Input file ATM_COR (see function ATMCOR(LAT) )			Page 32
Auxiliary programs									Page 33
	Program PP, version 4.00 (calculates production of calcite
		and Corg								Page 33
	Program REDUCE, version 2.00 (for data reduction)			Page 37
	Program TEST (for testing subroutines and functions)			Page 40




































	PROGRAM COCDIA
C
C	V. 10.04, created from COCDIA1003.for on 6.10.1994
C
C	One dimensional model to calculate growth of coccoliths
C	and diatoms. Includes production of attached coccoliths.
C	Detachment of liths is a function of the lith per cell ratio.
C	Includes free liths and the new light model for diatom growth.
C	Maximum growth rate of diatoms and coccolithophorids are 
C	calculated from Goldman and Carpenter, 1974.
C	Calculates annual carbonate production of E. huxleyi.
C	Includes expanded light model with light attenuation by
C	suspended liths.
C	Accelerated version of COCDIA310.
C	Includes routine for latitudinal mixed layer field.
C	Calculates annual organic carbon production.
C	Includes vertical and latitudinal variations of deep nitrate
C	and silicate fields.
C	Includes seasonal variation of mixing coefficients K1 and K2.
C	Calculates depth of photic zone. 
C	Preparation for output for annual production routines 
C	in external program PP200.
C	Includes seasonal variation of grazing and mortality rates.
C	Gives three output files PLANKTON.DAT, GROWTH.DAT and PHYSICS.DAT.
C	Uses instantaneous flux of solar radiation at sea surface instead
C	of daily averaged flux of solar variation.
C	Zooplankton grazing only during night.
C	Includes seasonal and latitudinal values for the sea-surface albedo.
C	Includes zonal atmospheric correction factors for light.
C	Includes a new temperature model for the growth rate of E. huxleyi:
C	  Diatoms and coccolithophorids follow different temperature
C	  characteristics.
C
C	Numerical solution by Euler method.
C
	INTEGER N,NMAX,NOUT
C
	DOUBLE PRECISION PM,PMN,DM,DMN,NM,NMN
	DOUBLE PRECISION SM,SMN
	DOUBLE PRECISION PT,PTN,DT,DTN,NT,NTN
	DOUBLE PRECISION ST,STN
	DOUBLE PRECISION PB,DB,NB,SB
C
	DOUBLE PRECISION FPM,FDM,FNM,FSM
	DOUBLE PRECISION FPT,FDT,FNT,FST
	DOUBLE PRECISION YMN,YTN
	DOUBLE PRECISION LENGTH,DELTAT,DLAT,TIME
C
	DOUBLE PRECISION HM,HT,HP,IM,IT
	DOUBLE PRECISION TM,TT,PHID,PHIC
	DOUBLE PRECISION GMAXE,GMAXD
	DOUBLE PRECISION ALPHAD,ALPHAC,CALC
C
	DOUBLE PRECISION LAM,LFM,LAMN,LFMN,FLAM,FLFM
	DOUBLE PRECISION LAT,LFT,LATN,LFTN,FLAT,FLFT
	DOUBLE PRECISION LAB,LFB
C
	DOUBLE PRECISION HMFELD(0:14),ALFELD(0:14)
	DOUBLE PRECISION A,B,C,D,DEPTH
C
	DOUBLE PRECISION ISURF,SLIGHT,LIGHTM,LIGHTT
	DOUBLE PRECISION TEMPM,TEMPT
	DOUBLE PRECISION MALPHAC,TALPHAC,MALPHAD,TALPHAD
	DOUBLE PRECISION PHICM,PHICT,PHIDM,PHIDT
	DOUBLE PRECISION GMXME,GMXMD,GMXTE,GMXTD
	DOUBLE PRECISION CALCM,CALCT
	DOUBLE PRECISION ATMCOR,ATRAN
C
	DOUBLE PRECISION KEXR,KEXG,KEXB
	DOUBLE PRECISION KEXRM,KEXGM,KEXBM
	DOUBLE PRECISION KEXRT,KEXGT,KEXBT
C
	CHARACTER*59 HEADER
C
C
C	Definition of parameters in include file COCDAT.inc:
C
	INCLUDE COCDAT.inc
	PARAMETER(NOUT=1)
C
C
	WRITE(9,*) '*** COCDIA Version 10.04 ***'
	WRITE(9,*) ' '
C
C	Read initial conditions from file "INITCON.DAT':
	OPEN (20,FILE='INITCON.DAT',STATUS='OLD')
C
	WRITE(9,*) 'Enter initial values for PM,LAM,LFM,DM,NM,SM 
     1  (=Mixed layer)'
	READ(20,*) PM,LAM,LFM,DM,NM,SM
	WRITE(9,*) PM,LAM,LFM,DM,NM,SM
	WRITE(9,*) 'Enter initial values for PT,LAT,LFT,DT,NT,ST
     2  (=Thermocline layer)'
     	READ(20,*) PT,LAT,LFT,DT,NT,ST
	WRITE(9,*) PT,LAT,LFT,DT,NT,ST
	WRITE(9,*) 'Enter values for PB,LAB,LFB,DB,NB,SB
     3  (=Bottom layer)'
	READ(20,*) PB,LAB,LFB,DB,NB,SB
	WRITE(9,*) PB,LAB,LFB,DB,NB,SB
C
C
C	Determine length of simulation (DELTAT=1/24)
C
	WRITE(9,*) 'The time-step DELTAT is 1 hour (i.e. 1/24)'
	WRITE(9,*) ' '
	WRITE(9,*) 'Take care for output of values:'
	WRITE(9,*) 'NOUT=1:  output every hour'
	WRITE(9,*) 'NOUT=12: output every 12 hours'
	WRITE(9,*) 'NOUT=24: output every 24 hours'
	WRITE(9,*) 'i.e. every day at midnight'
	WRITE(9,*) ' '	
	WRITE(9,*) 'Now, enter the length of the simulation'
	WRITE(9,*) '           (in days)'
C
	READ(9,*) LENGTH
C
	CLOSE (20)
C
C	Calculate maximum number of iterations (=NMAX)
C
	NMAX=INT(LENGTH/DELTAT)
C
	WRITE(9,*) 'Number of iterations= ',NMAX
C
	PAUSE 99
 99	CONTINUE
C
C	Initialize iteration (time increment is N).
C
	N=0
	TIME=N*DELTAT
C
C
C
C	Preparation to calculate mixed layer depths at latitude DLAT:
C	Reading Mixed-layer depth matrix (external file 'MIXED_LAYER_DEPTHS').
C	There are 7 latitudinal stations (running variable I in DO-loop 2).
C	The first column is the latitude, the next column is the starting
C	mixed layer (day 0), the next 12 columns are monthly average mixed layer
C	depths, the last column is the mixed layer depth at day 364).
C
	OPEN (16,FILE='MIXED_LAYER_DEPTHS',STATUS='OLD')
 250	READ(16,*,END=998) (HMFELD(I),I=0,14)
	IF (HMFELD(0).EQ.DLAT) THEN
	  CLOSE(16)
	  GOTO 100
	END IF
	GOTO 250
C
 100	CONTINUE
C
C
C 
C***	Now, preparing to calculate the sea-surface albedo
C	Reading Albedo data matrix and finding proper latitude:
C	In the file SEA_SURFACE_ALBEDO are the seasonal values
C	for the albedo at a given latitude (i.e. the albedo
C	at mid-time of each month). In the first column,
C	the latitude is stored.
C
	OPEN (16,FILE='SEA_SURFACE_ALBEDO',STATUS='OLD')
C
C	Reading header:
	READ(16,*) HEADER
C
 251	READ(16,*,END=997) (ALFELD(I),I=0,14)
	IF (ALFELD(0).EQ.DLAT) THEN
	  CLOSE (16)
	  GOTO 200
	END IF
	GOTO 251
C
 200	CONTINUE
C
C
C
C	Preparation to calculate the latitudinal reduction of light
C	during travelling through the atmosphere. ATRAN is the fraction
C	of solar radiation within PAR at the top of the atmosphere,
C	which eventually reaches the sea-surface (ATRAN stands for 
C	atmospheric transmission):
C
	ATRAN=ATMCOR(DLAT)
C
C
C
C	Preparation to calculate deep nutrient field at latitude DLAT.
C	Coefficients for nitrate are a and b, those for silicate are c and d:
C
	A=-7.0967D-4*DLAT+4.5158D-2
	B=0.54572D0*DLAT-17.366D0
	C=-3.1026D-4*DLAT+2.3291D-2
	D=0.30478D0*DLAT-10.9D0
C
C	Calculate initial deep nutrients at latitude DLAT for nitrate and silicate:
C
	DEPTH=HM(TIME,HMFELD)+HT
	NB=A*DEPTH+B
	SB=C*DEPTH+D
C
C
C
C	Now, calculate initial depth of the photic zone (HP):
C
	CALL PHOZM(PM,DM,LFM,HP)
C
C	In case that photic zone reaches deeper than the mixed layer:
C
	IF (HP.GT.HM(TIME,HMFELD)) THEN
	  CALL PHOZT(TIME,HMFELD,PM,DM,LFM,PT,DT,LFT,HP)
	END IF
C
C
C
C	Open three external files:
C	The file 'PLANKTON.DAT' contains standing stocks of plankton
C	populations and nutrient concentrations. The file 'PHYSICS.DAT.'
C	contains growth characteristics of diatoms and coccolithophorids,
C	and the file 'LIGHT.DAT' contains the depth of the mixed layer,
C	photic zone, light and temperature values.
C
	OPEN(15,FILE='PLANKTON.DAT',STATUS='NEW')
	OPEN(17,FILE='PHYSICS.DAT',STATUS='NEW')
	OPEN(18,FILE='LIGHT.DAT',STATUS='NEW')
C
C	Prepare for output of initial values by calculating common values:
C
	SLIGHT=ISURF(TIME,ALFELD,ATRAN)
	LIGHTM=IM(TIME,PM,DM,LFM,HMFELD,ALFELD,ATRAN)
	LIGHTT=IT(TIME,PT,DT,LFT,HMFELD,ALFELD,ATRAN)
	TEMPM=TM(TIME)
	TEMPT=TT(TIME)
	MALPHAC=ALPHAC(LIGHTM,TEMPM)
	TALPHAC=ALPHAC(LIGHTT,TEMPT)
	MALPHAD=ALPHAD(LIGHTM,TEMPM)
	TALPHAD=ALPHAD(LIGHTT,TEMPT)
	PHICM=PHIC(NM)
	PHICT=PHIC(NT)
	PHIDM=PHID(SM,NM)
	PHIDT=PHID(ST,NT)
	GMXME=GMAXE(TEMPM)
	GMXMD=GMAXD(TEMPM)
	GMXTE=GMAXE(TEMPT)
	GMXTD=GMAXD(TEMPT)
	KEXRM=KEXR(PM,DM)
	KEXGM=KEXG(PM,DM,LFM)
	KEXBM=KEXB(PM,DM,LFM)
	KEXRT=KEXR(PT,DT)
	KEXGT=KEXG(PT,DT,LFT)
	KEXBT=KEXB(PT,DT,LFT)
	CALCM=CALC(LIGHTM)
	CALCT=CALC(LIGHTT)
C
	WRITE(15,*) 'Latitude= ',DLAT,' DELTAT= ',DELTAT
	WRITE(15,*) 'Time (Days),PM,PT,DM,DT,LAM,LAT,LFM,LFT,NM,
     1NT,SM,ST,NB,SB'
	WRITE(15,*) TIME,PM,PT,DM,DT,LAM,LAT,LFM,LFT,NM,
     1NT,SM,ST,NB,SB
C
	WRITE(17,*) 'Latitude= ',DLAT,' DELTAT= ',DELTAT
	WRITE(17,*) 'Time(days),HM,HP,MALPHAC,TALPHAC,MALPHAD,
     1TALPHAD,PHICM,PHICT,PHIDM,PHIDT,GMXME,GMXMD,GMXTE,
     6GMXTD,CALCM,CALCT'
     	WRITE(17,*) TIME,HM(TIME,HMFELD),HP,MALPHAC,TALPHAC,
     2MALPHAD,TALPHAD,PHICM,PHICT,PHIDM,PHIDT,GMXME,GMXMD,
     3GMXTE,GMXTD,CALCM,CALCT
C
     	WRITE(18,*) 'TIME(days),ISURF,LIGHTM,LIGHTT,TEMPM,TEMPT,
     4KEXRM,KEXGM,KEXBM,KEXRT,KEXGT,KEXBT'
	WRITE(18,*) TIME,SLIGHT,LIGHTM,LIGHTT,TEMPM,TEMPT,KEXRM,
     5KEXGM,KEXBM,KEXRT,KEXGT,KEXBT
C
C
C	Iteration:
C	Calculation of components: The general form for each component is:
C
C	             XMN=XM+DELTAT*f(variables,N)
C
C	where f is a function, and N is time (Euler method). In order to
C	easily modify the program easily to a Runge-Kutta method, f is
C	calculated in different subroutines for each component (the 
C	subroutines return FXM for component X in the mixed layer and FXT
C	for component X in the thermocline layer).
C
C
	WRITE(9,*) '. . .Calculating. . .'
C
C	First calculate new deep nutrient concentration after time increment:
C
 5	TIME=N*DELTAT
 	DEPTH=HM(TIME,HMFELD)+HT
	NB=A*DEPTH+B
 	SB=C*DEPTH+D
C
C	Print message to screen every 24'th step while the program
C	is running:
C
 	IF (MOD(N,24).EQ.0.0) THEN
	  WRITE(9,*) '  Time: ',TIME
	END IF
C
C
C  Euler method:
C**** Mixed Layer:
	CALL MIXLAY(PM,PT,DM,DT,LAM,LAT,LFM,LFT,
     1	NM,NT,SM,ST,TIME,FPM,FDM,FNM,FSM,FLAM,FLFM,
     2  HMFELD,ALFELD,ATRAN)
C
	PMN=PM+DELTAT*FPM
	DMN=DM+DELTAT*FDM
	NMN=NM+DELTAT*FNM
	SMN=SM+DELTAT*FSM
	LAMN=LAM+DELTAT*FLAM
	LFMN=LFM+DELTAT*FLFM
C
C
C**** Thermocline-layer:
	CALL TERLAY(PM,PT,PB,DM,DT,DB,NM,NT,NB,
     1	SM,ST,SB,LAM,LAT,LAB,LFM,LFT,LFB,TIME,
     2	FPT,FDT,FNT,FST,FLAT,FLFT,HMFELD,
     3  ALFELD,ATRAN)
C
	PTN=PT+DELTAT*FPT
	DTN=DT+DELTAT*FDT
	NTN=NT+DELTAT*FNT
	STN=ST+DELTAT*FST
	LATN=LAT+DELTAT*FLAT
	LFTN=LFT+DELTAT*FLFT
C
C
C*****	Correction for entrainment or detrainment of material (cells, nutrients,
C	attached and free liths) as the mixed layer moves up or down. The new
C	values are stored in PMN,DMN,NMN,SMN,LAMN,LFMN and PTN,DTN,NTN,STN,LATN,
C	LFTN (the old values from above are simply written over):
C
C       Coccolithophorids:
	CALL ENTRAIN(TIME,PMN,PTN,PB,YMN,YTN,HMFELD)
	PMN=YMN
	PTN=YTN
C
C       Diatoms:
	CALL ENTRAIN(TIME,DMN,DTN,DB,YMN,YTN,HMFELD)
	DMN=YMN
	DTN=YTN
C
C       Nitrate:
	CALL ENTRAIN(TIME,NMN,NTN,NB,YMN,YTN,HMFELD)
	NMN=YMN
	NTN=YTN	
C
C       Silicate:
	CALL ENTRAIN(TIME,SMN,STN,SB,YMN,YTN,HMFELD)
	SMN=YMN
	STN=YTN
C
C	Attached liths:
	CALL ENTRAIN(TIME,LAMN,LATN,LAB,YMN,YTN,HMFELD)
	LAMN=YMN
	LATN=YTN
C
C	Free liths:
	CALL ENTRAIN(TIME,LFMN,LFTN,LFB,YMN,YTN,HMFELD)
	LFMN=YMN
	LFTN=YTN
C
C
C*****	Stop phytoplankton dying in winter:
	IF (PMN.LT.PB) THEN
	  PMN=PB
	ELSE IF (PTN.LT.PB) THEN
	  PTN=PB
	ELSE IF (DMN.LT.DB) THEN
	  DMN=DB
	ELSE IF (DTN.LT.DB) THEN
	  DTN=DB
	END IF
C
C****	Calculate depth of photic zone for determination of vertically integrated
C	annual production Corg and carbonate in external program PP. The depth
C	of the photic zone is calculated by the subroutines PHOZM and PHOZT, and is
C	stored in HP (in meters).
C
	CALL PHOZM(PM,DM,LFM,HP)
C
C	In case that photic zone reaches deeper than the mixed layer:
C
	IF (HP.GT.HM(TIME,HMFELD)) THEN
	  CALL PHOZT(TIME,HMFELD,PM,DM,LFM,PT,DT,LFT,HP)
	END IF
C
C
C
C**** Output of components to external files 'PLANKTON.DAT',
C	GROWTH.DAT and 'PHYSICS.DAT'.
C	Output to files every NOUT'th increment of N.
C
C	First calculate common seasonal variables:
C
	SLIGHT=ISURF(TIME,ALFELD,ATRAN)
	LIGHTM=IM(TIME,PM,DM,LFM,HMFELD,ALFELD,ATRAN)
	LIGHTT=IT(TIME,PT,DT,LFT,HMFELD,ALFELD,ATRAN)
	TEMPM=TM(TIME)
	TEMPT=TT(TIME)
	MALPHAC=ALPHAC(LIGHTM,TEMPM)
	TALPHAC=ALPHAC(LIGHTT,TEMPT)
	MALPHAD=ALPHAD(LIGHTM,TEMPM)
	TALPHAD=ALPHAD(LIGHTT,TEMPT)
	PHICM=PHIC(NM)
	PHICT=PHIC(NT)
	PHIDM=PHID(SM,NM)
	PHIDT=PHID(ST,NT)
	GMXME=GMAXE(TEMPM)
	GMXMD=GMAXD(TEMPM)
	GMXTE=GMAXE(TEMPT)
	GMXTD=GMAXD(TEMPT)
	KEXRM=KEXR(PM,DM)
	KEXGM=KEXG(PM,DM,LFM)
	KEXBM=KEXB(PM,DM,LFM)
	KEXRT=KEXR(PT,DT)
	KEXGT=KEXG(PT,DT,LFT)
	KEXBT=KEXB(PT,DT,LFT)
	CALCM=CALC(LIGHTM)
	CALCT=CALC(LIGHTT)
C
	IF (MOD(N,NOUT).EQ.0.0) THEN
C	++PLANKTON.DAT:
	  WRITE(15,*) (N+1)*DELTAT,PMN,PTN,DMN,DTN,LAMN,LATN,
     1LFMN,LFTN,NMN,NTN,SMN,STN,NB,SB
C
C	++PHYSICS.DAT:
     	WRITE(17,*) (N+1)*DELTAT,HM(TIME,HMFELD),HP,MALPHAC,
     3TALPHAC,MALPHAD,TALPHAD,PHICM,PHICT,PHIDM,PHIDT,
     4GMXME,GMXMD,GMXTE,GMXTD,CALCM,CALCT
C
C	++LIGHT.DAT:
	WRITE(18,*) (N+1)*DELTAT,SLIGHT,LIGHTM,LIGHTT,TEMPM,
     5TEMPT,KEXRM,KEXGM,KEXBM,KEXRT,KEXGT,KEXBT
C
	END IF
C
C
C	Write over old concentrations of components by the new values:
C
	PM=PMN
	DM=DMN
	NM=NMN
	SM=SMN
	PT=PTN
	DT=DTN
	NT=NTN
	ST=STN
	LAM=LAMN
	LAT=LATN
	LFM=LFMN
	LFT=LFTN
C
C	Stop criterium of iteration:
C
	IF (N+1.GE.NMAX) THEN
	  STOP
	END IF
C
C	Calculate the next time-step:
C
	N=N+1
	GOTO 5
C
 998	WRITE(9,*) 'End of MIXED_LAYER_DEPTH reached'
 	STOP
 997	WRITE(9,*) 'End of SEA_SURFACE_ALBEDO reached'
 	STOP
C
	END
	
	













	
	
C**************************************************************************
	SUBROUTINE MIXLAY(PM,PT,DM,DT,LAM,LAT,LFM,LFT,
     1	NM,NT,SM,ST,TIME,FPM,FDM,FNM,FSM,FLAM,FLFM,
     2  HMFELD,VALUE,TRANS)
C**************************************************************************
C
C	Calculates the components in the mixed layer due to production,
C	uptake, loss and mixing.
C
C	Input values into the subroutine are:
C	PM,PT,DM,DT:		mg Chla m-3
C	LAM,LAT,LFM,LFT:	LITH m-3
C	NM,NT,SM,ST:		mM m-3
C	TIME:			days
C
C	Output values are:
C	FPM,FDM:		mg Chla m-3 day-1
C	FLAM,FLFM:		Lith m-3 day-1
C	FNM,FSM:		mM m-3 day-1
C
	DOUBLE PRECISION PM,PT,DM,DT,LAM,LAT,LFM,LFT
	DOUBLE PRECISION NM,NT,SM,ST,TIME
	DOUBLE PRECISION FPM,FDM,FNM,FSM,FLAM,FLFM
	DOUBLE PRECISION HM,IM,TM,RDET,ALPHAC,ALPHAD
	DOUBLE PRECISION GMAXD,GMAXE,MAXGD,MAXGE
	DOUBLE PRECISION PHIC,PHID,CALC
	DOUBLE PRECISION DEPTH,LIGHT,TEMP,DETACH
	DOUBLE PRECISION GROWC,GROWD
	DOUBLE PRECISION MCM,MDM,GRAZM,GRAZD
	DOUBLE PRECISION VCOC,VDIA,GAMMAC
	DOUBLE PRECISION GAMMAD,EPSCOC,EPSDIA,CHLPC,K2,K2S
	DOUBLE PRECISION HMFELD(0:14),VALUE(0:14),TRANS
C
	INCLUDE COCDAT.inc
C
C	Calculate seasonal variations:
C
	DEPTH=HM(TIME,HMFELD)
	LIGHT=IM(TIME,PM,DM,LFM,HMFELD,VALUE,TRANS)
	TEMP=TM(TIME)
	K2=K2S*(DEPTH/25.0D0)**1.5
	MAXGD=GMAXD(TEMP)
	MAXGE=GMAXE(TEMP)
	DETACH=RDET(PM,LAM)
	GROWC=ALPHAC(LIGHT,TEMP)*MAXGE*PHIC(NM)
	GROWD=ALPHAD(LIGHT,TEMP)*MAXGD*PHID(SM,NM)
C
C
C	Calculate grazing rates in a day-night cycle:
C	No grazing in mixed layer during day (grazers are at great depths),
C	and only during night time grazers are active.
C
	IF (LIGHT.EQ.0.0) THEN
	  GRAZC=MCM(TIME)
	  GRAZD=MDM(TIME)
	ELSE
	  GRAZC=0.0
	  GRAZD=0.0
	END IF
C
C
C	Now, calculate individual constituents:
C	Coccolithophorids:
C
	FPM=(GROWC-GRAZC-VCOC/DEPTH)*PM+
     1	K2*(PT-PM)/DEPTH
C
C	Diatoms:
C
	FDM=(GROWD-GRAZD-VDIA/DEPTH)*DM+
     1	K2*(DT-DM)/DEPTH
C
C	Nitrate:
C
	FNM=-GAMMAC*(GROWC-EPSCOC*GRAZC)*PM-
     1	GAMMAD*(GROWD-EPSDIA*GRAZD)*DM+
     2  K2*(NT-NM)/DEPTH
C
C	Silicate:
C
	FSM=-GAMMAD*GROWD*DM+K2*(ST-SM)/DEPTH
C
C	Attached liths:
C
	FLAM=(CALC(LIGHT)-DETACH)*CHLPC*PM-
     1  GRAZC*LAM-VCOC*LAM/DEPTH+K2*(LAT-LAM)/DEPTH
C
C	Free liths:
C
	FLFM=DETACH*CHLPC*PM-GRAZC*LFM+
     1	K2*(LFT-LFM)/DEPTH
C
	RETURN
	END







C**************************************************************************
	SUBROUTINE TERLAY(PM,PT,PB,DM,DT,DB,NM,NT,NB,
     1	SM,ST,SB,LAM,LAT,LAB,LFM,LFT,LFB,TIME,
     2	FPT,FDT,FNT,FST,FLAT,FLFT,HMFELD,VALUE,TRANS)
C**************************************************************************
C
C	Calculates the components in the thermocline layer due to production,
C	uptake, loss and mixing.
C
C	Input values into the subroutine are:
C	PM,PT,DM,DT:		mg Chla m-3
C	LAM,LAT,LFM,LFT:	LITH m-3
C	NM,NT,SM,ST:		mM m-3
C	TIME:			days
C
C	Output values are:
C	FPM,FDM:		mg Chla m-3 day-1
C	FLAM,FLFM:		Lith m-3 day-1
C	FNM,FSM:		mM m-3 day-1
C
	DOUBLE PRECISION PM,PT,PB,DM,DT,DB,NM,NT,NB
	DOUBLE PRECISION SM,ST,SB,LAM,LAT,LAB
	DOUBLE PRECISION LFM,LFT,LFB,TIME
	DOUBLE PRECISION FPT,FDT,FNT,FST,FLAT,FLFT
	DOUBLE PRECISION IT,TT,RDET,ALPHAC,ALPHAD
	DOUBLE PRECISION GMAXD,GMAXE,MAXGD,MAXGE
	DOUBLE PRECISION CALC,PHIC,PHID
	DOUBLE PRECISION LIGHT,TEMP,DETACH,GROWC
	DOUBLE PRECISION GROWD,VCOC,VDIA
	DOUBLE PRECISION MCT,MDT,GRAZC,GRAZD
	DOUBLE PRECISION GAMMAC,GAMMAD,EPSCOC,EPSDIA
	DOUBLE PRECISION K1,K2,K1S,K2S,HT,CHLPC
	DOUBLE PRECISION HMFELD(0:14),VALUE(0:14),HM
	DOUBLE PRECISION DEPTH,TRANS
C
	INCLUDE COCDAT.inc
C
C	Calculate seasonal variations:
C
	DEPTH=HM(TIME,HMFELD)
	LIGHT=IT(TIME,PT,DT,LFT,HMFELD,VALUE,TRANS)
	TEMP=TT(TIME)
	K2=K2S*(DEPTH/25.0D0)**1.5
	K1=K1S*(DEPTH/25.0D0)**1.5
	MAXGD=GMAXD(TEMP)
	MAXGE=GMAXE(TEMP)
	DETACH=RDET(PT,LAT)
	GROWC=ALPHAC(LIGHT,TEMP)*MAXGE*PHIC(NT)
	GROWD=ALPHAD(LIGHT,TEMP)*MAXGD*PHID(ST,NT)
C
C
C	Simulate grazing rates in a day-night cycle:
C	No grazing in mixed layer during day (grazers are at great depths),
C	and only during night time grazers are active.
C
	IF (LIGHT.EQ.0.0) THEN
	  GRAZC=MCT(TIME)
	  GRAZD=MDT(TIME)
	ELSE
	  GRAZC=0.0
	  GRAZD=0.0
	END IF
C
C
C	Now, calculate individual constituents:
C	Coccolithophorids:
C
	FPT=(GROWC-GRAZC)*PT-
     1	(VCOC*(PT-PM)-K1*(PB-PT)+K2*(PT-PM))/HT
C
C	Diatoms:
C
	FDT=(GROWD-GRAZD)*DT-
     1	(VDIA*(DT-DM)-K1*(DB-DT)+K2*(DT-DM))/HT
C
C	Nitrate:
C
	FNT=-GAMMAC*(GROWC-EPSCOC*GRAZC)*PT-
     1	GAMMAD*(GROWD-EPSDIA*GRAZD)*DT+
     2  (K1*(NB-NT)-K2*(NT-NM))/HT
C
C	Silicate:
C
	FST=-GAMMAD*GROWD*DT+
     1	(K1*(SB-ST)-K2*(ST-SM))/HT
C
C	Attached liths:
C
	FLAT=(CALC(LIGHT)-DETACH)*CHLPC*PT-
     1	GRAZC*LAT-(VCOC*(LAT-LAM)-
     2	K1*(LAB-LAT)+K2*(LAT-LAM))/HT
C
C	Free liths:
C
	FLFT=DETACH*PT*CHLPC-GRAZC*LFT+
     1	(K1*(LFB-LFT)-K2*(LFT-LFM))/HT
C
	RETURN
	END
	

	







C**************************************************************************
	DOUBLE PRECISION FUNCTION ALPHAC(LIGHT,TEMP)
C**************************************************************************
C
C	Calculates the growth of coccolithophorids as a function of light 
C	according to Balch et al. (1992). Light enters the function in Wm-2.
C
C	PCMAX is in fmol C cell-1 hr-1 at 15 degrees celcius and determines
C	the efficiecy of C uptake (from Balch et al., 1992).
C	ALPHAC approaches to 1 under light saturation and has no units.
C
	DOUBLE PRECISION ALPHAC,LIGHT,I
	DOUBLE PRECISION GMAXE,TEMP
	DOUBLE PRECISION AP,CCPC
C
	INCLUDE COCDAT.inc
C
C	Conversion of Wm-2 into Ein m-2 s-1 (at 550nm):
C
	I=LIGHT/2.17541D5
C
	ALPHAC=1.0-EXP(-AP*3600.0*24.0*I/(CCPC*GMAXE(TEMP)))
C
	RETURN
	END



	
C**************************************************************************
	DOUBLE PRECISION FUNCTION ALPHAD(LIGHT,TEMP)
C**************************************************************************
C
C	Calculates the light dependent growth of diatoms as a function
C	of light. Light dependency in form of Platt et al., 1985.
C	Light enters the function in Wm-2 and is converted into
C	Ein m-2 s-1 for the calculation. ALPHAD approaches to 1 under light
C	saturation and has no units.
C
	DOUBLE PRECISION ALPHAD,LIGHT,I
	DOUBLE PRECISION GMAXD,TEMP
	DOUBLE PRECISION APD,CCPCD
C
	INCLUDE COCDAT.inc
C
C	Conversion of Wm-2 into Ein m-2 s-1 (at 550nm):
C
	I=LIGHT/2.17541D5
C
C	Now calculation of the specific growth rate:
C
	ALPHAD=1.0-EXP(-APD*3600.0*24.0*I/(CCPCD*GMAXD(TEMP)))
C
	RETURN
	END
	














C**************************************************************************
	DOUBLE PRECISION FUNCTION GMAXD(TEMP)
C**************************************************************************
C
C	Model G:
C	Calculates the maximum possible growth rate of diatoms according
C	to the Model of Goldman and Carpenter, 1974. TEMP enters the
C	function in degrees celsius. GMAXD is returned to the calling unit
C	in day-1 (=specific growth rate).
C
	DOUBLE PRECISION GMAXD,TEMP
C
	GMAXD=5.35D9*EXP(-6472.0/(TEMP+273.0))
C
	RETURN
	END
	
	
		
	

C**************************************************************************
	DOUBLE PRECISION FUNCTION GMAXE(TEMP)
C**************************************************************************
C
C	Subroutine to calculate the maximum specific growth rate of
C	E. huxleyi as a function of temperature. Data were from Dave Lesley,
C	Roger Harris and Maureen Conte (Poster during fifth GEM Meeting at
C	Blagnac, September 1994). TEMP enters the subroutine in degrees
C	Celsius, and GMAXE is returned to the calling unit in day-1.
C	GMAXE stands for GMAX for E miliania huxleyi.
C
	INTEGER I,K,SI,EI
	DOUBLE PRECISION T(1:10),G(1:10)
	DOUBLE PRECISION GMAXE,TEMP,SLOPE
C
C	Reading temperatures of curve:
C
	DATA (T(I),I=1,10) /5.0D0,6.0D0,9.0D0,12.0D0,
     *  15.0D0,18.0D0,21.0D0,24.0D0,27.0D0,30.0D0/
C
C	Reading associated growth rates of curve:
C
	DATA (G(I),I=1,10) /0.0D0,0.144D0,0.243D0,
     *  0.466D0,0.592D0,0.773D0,1.024D0,0.981D0,0.935D0,0.0D0/
C
C	Calculating specific growth rates:
C	First, growth rates outside temperature tolerance range
C	for E. huxleyi are zero:
C
	IF ((TEMP.LT.T(1)).OR.(TEMP.GT.T(10))) THEN
	  GMAXE=0.0D0
	  RETURN
	END IF
C
C	Now, calculating growth rates if actual temperature is inside the
C	temperature tolerance range for E. huxleyi as function of TEMP:
C
	  DO 10, K=2,10
	    IF (TEMP.LE.T(K)) THEN
	      SI=K-1
	      EI=K
	      SLOPE=(G(EI)-G(SI))/(T(EI)-T(SI))
	      GMAXE=G(SI)+(TEMP-T(SI))*SLOPE
	      RETURN
	    END IF
 10	  CONTINUE
C
	END

C**************************************************************************
	DOUBLE PRECISION FUNCTION PHIC(NITRATE)
C**************************************************************************
C
C	Calculates the nitrate dependent growth for coccolithophorids
C	according to Michaelis Menten. The values of PHIC range from
C	0 to 1 and have no dimensions.
C
	DOUBLE PRECISION PHIC,NITRATE
	DOUBLE PRECISION NCHALF
	DOUBLE PRECISION VAR1
C
	INCLUDE COCDAT.inc
C
	VAR1=NITRATE+NCHALF
	PHIC=NITRATE/VAR1
C
	RETURN
	END
	
		
C**************************************************************************
	DOUBLE PRECISION FUNCTION PHID(SILICATE,NITRATE)
C**************************************************************************
C
C	Calculates the nitrate and silicate dependent growth term for
C	diatoms according to Liebig's law of the minimum.
C
	DOUBLE PRECISION PHID,SILICATE,NITRATE
	DOUBLE PRECISION PHINIT,PHISIL,NDHALF,SHALF
C
	INCLUDE COCDAT.inc
C
	PHINIT=NITRATE/(NITRATE+NDHALF)
	PHISIL=SILICATE/(SILICATE+SHALF)
C
C	Calculate the smaller of the two:
C
	IF (PHINIT.LE.PHISIL) THEN
	  PHID=PHINIT
	ELSE IF (PHINIT.GT.PHISIL) THEN
	  PHID=PHISIL
	END IF
C
	RETURN
	END
	
	
	
C*************************************************************************
	DOUBLE PRECISION FUNCTION CALC(LIGHT)
C*************************************************************************
C
C	Calculates the production of coccoliths per cell as a function
C	of light, according to Balch et al., 1992. LIGHT enters the
C	function in Wm-2, and CALC is returned to the calling unit in
C	Lith cell-1 day-1. CCMAX is the maximum calcite carbon uptake
C	rate in fmolC cell-1 hr-1 and CDAR is the dark calcification rate
C	also in fmolC cell-1 hr-1.
C
	DOUBLE PRECISION CALC,LIGHT
	DOUBLE PRECISION CCMAX,RACC,CON,CDAR
C
	PARAMETER(CON=6.942554D-1,RACC=7.0682D-3)
	PARAMETER(CCMAX=26.0D0,CDAR=0.0D0)
C
	CALC=CON*(CCMAX*(1.0-EXP(-RACC*LIGHT))+CDAR)
C
	RETURN
	END
C*************************************************************************
	DOUBLE PRECISION FUNCTION RDET(EHUX,LITH)
C*************************************************************************
C
C	Calculates the rate of detachment of liths as a function of the
C	ratio of attached liths per cell. EHUX enters the function in
C	mgChla m-3, and LITH in Lith m-3. RDET is returned to the calling
C	unit in Lith cell-1 day-1.
C
	DOUBLE PRECISION RDET,EHUX,LITH
	DOUBLE PRECISION CHLPC,OMEGA,QMAX,PI
	DOUBLE PRECISION Q
C
C	Omega is a constant and indicates how fast RDET reaches a maximum
C	value. Q is the actual ratio of liths per cell in Lith cell-1 and
C	QMAX is a maximal ratio of liths per cell a coccosphere may reach, 
C	and is of the order of 80 Liths cell-1 (Balch et al., 1993).
C
	INCLUDE COCDAT.inc
	PARAMETER(OMEGA=1.0D0,QMAX=80.0D0)
C
	Q=LITH/(EHUX*CHLPC)
	RDET=OMEGA*TAN(PI*Q/(2.0*QMAX))
C
	RETURN
	END





C*************************************************************************
	DOUBLE PRECISION FUNCTION HM(TIME,HEIGHT)
C*************************************************************************
C
C	Calculates the thickness of the mixed layer in meters as a 
C	function of time (in days). The data use seasonal values for the top
C	of the thermocline, which were read from maps given in Robinson, Bauer
C	and Schroeder, 1979. The necessary data are stored in the external
C	file MIXED_LAYER_DEPTH.DAT. A modified method from Mike Ainsworth,
C	PML, was used to calculate mixed layer depths in between the points.
C
	INTEGER I,K,SI,EI
	DOUBLE PRECISION DAY(1:14),HEIGHT(0:14)
	DOUBLE PRECISION HM,TIME,TR,SLOPE
C
	INCLUDE COCDAT.inc
C
C	In DAY(I) are the mid-month day numbers stored. In DAY(1) is the
C	first January (=day number 0), and in DAY(14) is the 31st December
C	(=day number 364) stored.
C
	DATA (DAY(I),I=1,14)/0.0D0,15.0D0,45.0D0,74.0D0,
     1	104.0D0,136.0D0,165.0D0,196.0D0,227.0D0,
     2  257.0D0,288.0D0,318.0D0,349.0D0,364.0D0/
C
C	First reduce multi-year cycles to a one-year cycle:
C	TR means the "time reduced". This implies, that there are
C	no inter-annual changes in the course of the mixed-layer depth.
C
	TR=MOD(TIME,364.0D0)
C
C	Now, find out to which interval the resetted time belongs:
C	(SI means Start-interval, EI means End-interval)
C
	DO 10, K=2,15
	  IF (TR.LE.DAY(K)) THEN
	    SI=K-1
	    EI=K
	    SLOPE=(HEIGHT(EI)-HEIGHT(SI))/(DAY(EI)-DAY(SI))
	    HM=HEIGHT(SI)+(TR-DAY(SI))*SLOPE
C
	    RETURN
	  END IF
 10	CONTINUE
C
	END






C*************************************************************************
	DOUBLE PRECISION FUNCTION TM(TIME)
C*************************************************************************
C
C	Calculates the temperature of the mixed layer as a function of
C	time.
C
	DOUBLE PRECISION TM,TIME,TMAXM,TMINM
	DOUBLE PRECISION PI,ARG
C
	INCLUDE COCDAT.inc
C
	ARG=2.0*PI*(TIME-135.75)/365.0
	TM=TMINM+0.5*(TMAXM-TMINM)*
     3  (1.0+SIN(ARG))
C
	RETURN
	END
	
	
	
C*************************************************************************
	DOUBLE PRECISION FUNCTION TT(TIME)
C*************************************************************************
C
C	Calculates the temperature of the thermocline layer as a function
C	of time.
C
	DOUBLE PRECISION TT,TM,TB,TIME
C
	INCLUDE COCDAT.inc
C
	TT=0.5*(TM(TIME)+TB)
C
	RETURN
	END
	
	
	

C*************************************************************************
	DOUBLE PRECISION FUNCTION KEXR(PX,DX)
C*************************************************************************
C
C	Calculates the extinction koefficient for red light as a function
C	of chlorophyll concentration of coccolithophorids (PX) and diatoms
C	(DX). PX and DX enters in mgChla m-3, and KEXR is returned to the
C	calling unit in m-1.
C
	DOUBLE PRECISION KEXR,PX,DX,KR0
C
	INCLUDE COCDAT.inc
C
	KEXR=KR0+0.016*(PX+DX)
C
	RETURN
	END	
C*************************************************************************
	DOUBLE PRECISION FUNCTION KEXG(PX,DX,LFX)
C*************************************************************************
C
C	Calculates the extinction koefficient for green light as a function
C	of chlorophyll concentration of coccolithophords (PX), diatoms
C	(DX) and of isolate liths (LFX). PX and DX enter in mg Chla m-3,
C	LFX enters in Lith m-3, and KEXG returns to the calling unit in m-1.
C	The C/Chla=40 in diatoms. Backscatter of light due to coccoliths 
C	according to Balch et al., 1991.
C
	DOUBLE PRECISION KEXG,PX,DX,LFX
	DOUBLE PRECISION CHLPC,CCPCD
C
	INCLUDE COCDAT.inc
C
	KEXG=7.028D-2+1.71D-11*(PX*CHLPC+DX/(CCPCD*3.00287D-13))+
     1  1.2896D-13*LFX-8.4063D-26*LFX**2
C
	RETURN
	END
	
	
	
	
C*************************************************************************
	DOUBLE PRECISION FUNCTION KEXB(PX,DX,LFX)
C*************************************************************************
C
C	Calculates the extinction koefficient for blue light as a function
C	of chlorophyll concentration of coccolithophorids (PX) and diatoms
C	(DX) and of isolate liths (LFX). PX and DX enter in mg Chla m-3,
C	LFX enters in Lith m-3, and KEXB returns to the calling unit in m-1.
C	The C/Chla=40 in diatoms. Backscatter of light due to coccoliths 
C	according to Balch et al., 1991.
C
	DOUBLE PRECISION KEXB,PX,DX,LFX
	DOUBLE PRECISION CHLPC,CCPCD
C
	INCLUDE COCDAT.inc
C
C	Calculation of chlorophyll content per cell in diatoms:
C
	KEXB=2.7836D-2+5.869D-11*(PX*CHLPC+
     1  DX/(CCPCD*3.00287D-13))+
     2  1.4068D-13*LFX-5.2696D-26*LFX**2
C
	RETURN
	END





C*************************************************************************
	DOUBLE PRECISION FUNCTION IM(TIME,PX,DX,LFX,
     1  HMFELD,VALUE,TRANS)
C*************************************************************************
C
C	Calculates the average light intensity in the mixed layer with
C	thickness HM, and as a function of coccolithophorid- and diatom
C	chlorophyll concentrations PX and DX, and of concentrations of
C	isolated liths LFX, respectively. IM returns to the calling unit
C	in Wm-2.
C
	DOUBLE PRECISION IM,TIME,PX,DX,LFX
	DOUBLE PRECISION HM,ISURF,KEXR,KEXG,KEXB
	DOUBLE PRECISION VAR1,VAR2,VAR3
	DOUBLE PRECISION DEPTH,RED,GREEN,BLUE
	DOUBLE PRECISION HMFELD(0:14),VALUE(0:14),TRANS
C
	DEPTH=HM(TIME,HMFELD)
	RED=KEXR(PX,DX)
	GREEN=KEXG(PX,DX,LFX)
	BLUE=KEXB(PX,DX,LFX)
C
	VAR1=(1-EXP(-RED*DEPTH))/RED
	VAR2=(1-EXP(-GREEN*DEPTH))/GREEN
	VAR3=(1-EXP(-BLUE*DEPTH))/BLUE
C
C	The entering surface light (ISURF) here is PAR only. 
C	(The 50% loss of irradiation, which is infrared, has already been
C	accounted for in function ISURF).
C
	IM=ISURF(TIME,VALUE,TRANS)*
     1  (VAR1+VAR2+VAR3)/(3.0*DEPTH)
C
	RETURN
	END
	



C*************************************************************************
	DOUBLE PRECISION FUNCTION IT(TIME,PX,DX,LFX,
     1  HMFELD,VALUE,TRANS)
C*************************************************************************
C
C	Calculates the average light intensity in the thermocline layer in
C	the presence of coccoliths (PX), diatoms (DX) and isolate liths (LFX).
C	PX and DX enter in mg Chla m-3, LFX in Lith m-3. IT returns to the
C	calling unit in Wm-2.
C
	DOUBLE PRECISION IT,TIME,PX,DX,LFX
	DOUBLE PRECISION HM,KEXR,KEXG,KEXB,ISURF
	DOUBLE PRECISION RED,GREEN,BLUE,DEPTH
	DOUBLE PRECISION VAR1,VAR2,VAR3,HT
	DOUBLE PRECISION HMFELD(0:14),VALUE(0:14)
	DOUBLE PRECISION TRANS
C
	INCLUDE COCDAT.inc
C
	DEPTH=HM(TIME,HMFELD)
	RED=KEXR(PX,DX)
	GREEN=KEXG(PX,DX,LFX)
	BLUE=KEXB(PX,DX,LFX)
C
	VAR1=(1.0-EXP(-HT*RED))*EXP(-DEPTH*RED)/RED
	VAR2=(1.0-EXP(-HT*GREEN))*EXP(-DEPTH*GREEN)/GREEN
	VAR3=(1.0-EXP(-HT*BLUE))*EXP(-DEPTH*BLUE)/BLUE
C
C	The entering surface light (ISURF) here is PAR only. 
C	(The 50% loss of irradiation, which is infrared, has already been
C	accounted for in function ISURF).
C
	IT=ISURF(TIME,VALUE,TRANS)*
     1  (VAR1+VAR2+VAR3)/(3.0*HT)
C
	RETURN
	END

	






	
	
C***************************************************	
	SUBROUTINE ENTRAIN(TIME,XM,XT,XB,XMN,XTN,HMFELD)
C*************************************************************************
C
C	Created 21.11.93, modified 24.11.93
C	Entrainment and detrainment in the mixed and thermocline layers:
C	The subroutine calculates the new concentration of component X 
C	after the mixed layer has moved up or down. M relates to the mixed
C	layer, T to the thermocline layer and B to the bottom layer. The
C	new value is returned to the calling unit in XMN.
C
	DOUBLE PRECISION TIME,DELTAT,XM,XT,XB,XMN,XTN
	DOUBLE PRECISION HM,HT,HMFELD(0:14)
	DOUBLE PRECISION DELHM,VAR1,DEPB,DEPN
C
	INCLUDE COCDAT.inc
C
	DEPB=HM(TIME,HMFELD)
	DEPN=HM(TIME+DELTAT,HMFELD)
C
	VAR1=DEPB/DEPN
	DELHM=DEPN-DEPB
C
C	Thermocline moving downward:
	IF (DELHM.GE.0.0D0) THEN
	  XMN=VAR1*XM+(DELHM*XT)/DEPN
	  XTN=XT-DELHM*(XT-XB)/HT
C	Thermocline moving upward:
	ELSE IF (DELHM.LT.0.0D0) THEN
	  XMN=XM
	  XTN=XT-DELHM*(XM-XT)/HT
	END IF	
C
	RETURN
	END
	
	
	



C*************************************************************************	
	DOUBLE PRECISION FUNCTION ISURF(TIME,VALUE,TRANS)
C*************************************************************************
C
C	Calculates the instantaneous flux of PAR of the solar radiation,
C	which penetrates the sea surface during the year and at a given
C	latitude. TIME enters the function in days and the light intensity
C	is returned to the calling unit in Wm-2. This function needs a time
C	increment of DELTAT=1/24 (i.e. one hour). After Sellers, 1965.
C
	INTEGER N
	DOUBLE PRECISION ISURF,TIME,LAT,DLAT,DELTAT
	DOUBLE PRECISION HALDAY,DIST,DAV,S,DEKL
	DOUBLE PRECISION PI,HOURCL,HOURRD,H
	DOUBLE PRECISION VAR1,VAR2,VAR3
	DOUBLE PRECISION VAR4,ALBEDO,VALUE(0:14),TRANS
C
	PARAMETER(DAV=1.496D11,S=1373.0D0)
C
C	S is the solar constant outside of the earth atmosphere in Wm-2
C	(after Kirk, 1994)
C
	INCLUDE COCDAT.inc
C
C	Convert DLAT into radians:
C
	LAT=PI*DLAT/180.0D0
C
C	Calculate hour during every day:
C
	N=INT(TIME/DELTAT)
C
C	HOURCL:= Time, in clock hours
C	HOURRD:= Time, in radians
C
	HOURCL=MOD(N,24)
	HOURRD=PI*(-1.0+HOURCL/12.0)
C
C	Calculate half-day length H (in radians):
C
	H=HALDAY(TIME)
C
C	Calculate instantaneous solar flux:
C
	IF (ABS(HOURRD).LT.ABS(H)) THEN
C
C	  Solar flux during day-time:
C
	  VAR1=SIN(LAT)*SIN(DEKL(TIME))
	  VAR2=COS(LAT)*COS(DEKL(TIME))*COS(HOURRD)
	  VAR3=S*(DAV/DIST(TIME))**2
	  VAR4=1.0-ALBEDO(TIME,VALUE)
C
	  ISURF=0.38*TRANS*VAR3*VAR4*(VAR1+VAR2)
C
C	  The constant factor of 0.38 accounts for the fraction of PAR (400nm-700nm)
C	  in the extraterrestrial radiation (38 %).
C	  In TRANS is stored the latitudinal relative fraction of light (with respect
C	  to the solar radiation at the top of the atmosphere), which reaches the
C	  sea-surface after passing the atmosphere.
C	  ALBEDO contains the seasonal and latitudinal value for the sea-surface albedo.
C
	ELSE
C
C	  Solar flux during night-time:
C
	  ISURF=0.0D0
	END IF
C
	RETURN
	END





C*************************************************************************
	DOUBLE PRECISION FUNCTION ALBEDO(TIME,VALUE)
C*************************************************************************
C
C	Calculates the variation of the sea-surface albedo as a 
c	function of time (in days) at a given latitude. The data
C	to approximate the seasonal curve are from Goldsmith and Bunker, 1979.
C	TIME enters in days, VALUE in %, and the value of ALBEDO is returned
C	to the calling unit in per cent. Method adapted from Mike Ainsworth,
C	PML.
C
	INTEGER I,K,SI,EI
	DOUBLE PRECISION DAY(1:14),VALUE(0:14)
	DOUBLE PRECISION ALBEDO,TIME,TR,SLOPE
C
	INCLUDE COCDAT.inc
C
C	The mid-month day numbers are stored in DAY(I). The value of
C	DAY(1) and DAY(14) are the first day in January, and the last
C	day during December, respectively.
C
	DATA (DAY(I),I=1,14)/0.0D0,15.0D0,45.0D0,74.0D0,
     1	104.0D0,136.0D0,165.0D0,196.0D0,227.0D0,
     2  257.0D0,288.0D0,318.0D0,349.0D0,364.0D0/
C
C	First reduce multi-year cycles to a single-year cycle (periodicity
C	of 1 year): TR means the "time reduced". This implies, that there are
C	no inter-annual changes in the course of the mixed-layer depth.
C
	TR=MOD(TIME,364.0D0)
C
C	Now, find out to which interval the resetted time belongs:
C	(SI means Start-interval, EI means End-interval). VALUE(SI) 
C	and VALUE(EI) contain the albedo at the begin or the end of
C	a particular interval between the mid-month dates.
C
	DO 10, K=2,15
	  IF (TR.LE.DAY(K)) THEN
	    SI=K-1
	    EI=K
	    SLOPE=(VALUE(EI)-VALUE(SI))/(DAY(EI)-DAY(SI))
	    ALBEDO=VALUE(SI)+(TR-DAY(SI))*SLOPE
C
	    RETURN
	  END IF
 10	CONTINUE
C
	END






C*************************************************************************	
	DOUBLE PRECISION FUNCTION HALDAY(TIME)
C*************************************************************************
C
C	Calculates the half day length at a location with
C	latitude L and at time N*DELTAT (in days, N=0 at January 1.
C	TIME enters the function in days, the result is returned to
C	the calling unit in radians.
C	After Sellers, 1965.
C
	DOUBLE PRECISION HALDAY,DEKL,LAT,DLAT,TIME
	DOUBLE PRECISION ARG,PI
C
	INCLUDE COCDAT.inc
C
C	Convert DLAT into radians:
	LAT=PI*DLAT/180.0D0

	ARG=-TAN(LAT)*TAN(DEKL(TIME))
C
	IF ((ARG).GE.1.0) THEN
	  HALDAY=0.0
	  RETURN
	ELSE IF ((ARG).LE.-1.0) THEN
	  HALDAY=PI
	  RETURN
	END IF
C
	HALDAY=DACOS(ARG)
C
	RETURN
	END
	
	





C*************************************************************************	
	DOUBLE PRECISION FUNCTION DEKL(TIME)
C*************************************************************************
C
C	Function to calculate the declination of the sun as a function
C	of time (N*DELTAT, in days). N is 0 on January 1 and 364 on December 31.
C	DEKL returns the declination in radians (im Bogenmass).
C
	DOUBLE PRECISION DEKL,TIME,PI
	DOUBLE PRECISION ARG
C
	INCLUDE COCDAT.inc
C
	ARG=2.0D0*PI*TIME/364.0D0-1.3831615D0
	DEKL=23.45D0*PI*SIN(ARG)/180.0D0
C
	RETURN
	END






C*************************************************************************	
	DOUBLE PRECISION FUNCTION DIST(TIME)
C*************************************************************************
C
C	Function to calculate the instantanous distance of the earth from
C	the sun as a function from time (T, in days). T is 0 on January 1
C	and 364 on December 31. The value of DIST is returned to the
C	calling unit in meters. From Mike Ainsworth, PML.
C
	DOUBLE PRECISION DIST,TIME,PI,DMIN,DMAX
	DOUBLE PRECISION ARG
C
	PARAMETER(DMIN=1.471D11,DMAX=1.521D11)
C	DMIN is the minimum distance of the earth from the sun in m,
C	DMAX is the maximum distance of the earth from the sun in m.
C
	INCLUDE COCDAT.inc
C
	ARG=(TIME-94.75)/365.0*2.0*PI
	DIST=DMIN+(DMAX-DMIN)/2.0*(SIN(ARG)+1.0)
C
	RETURN
	END





	
	
C*************************************************************************
	SUBROUTINE PHOZM(PX,DX,LFX,DEPTH)
C*************************************************************************
C
C	This subroutine calculates the depth of the photic zone inside of
C	the mixed layer. Intervals of calculation are 1 meter. The base of the
C	photic zone is considered at that depth, where 1% of the surface PAR
C	is attained. LPZ is the relative light intensity at the base of the
C	photic zone. PX,DX enter in units of mg Chla m-3, LFX in Lith m-3, and
C	Depth is returned to the calling unit in meters.
C
	DOUBLE PRECISION DEPTH,LPZ,PX,DX,LFX
	DOUBLE PRECISION KEXR,KEXG,KEXB
	DOUBLE PRECISION VAR1,VAR2,VAR3
C
	DEPTH=0.0D0
 5	VAR1=EXP(-DEPTH*KEXR(PX,DX))
	VAR2=EXP(-DEPTH*KEXG(PX,DX,LFX))
	VAR3=EXP(-DEPTH*KEXB(PX,DX,LFX))
	LPZ=VAR1+VAR2+VAR3
C
	IF (LPZ.LE.0.03) THEN
	  RETURN
	END IF
C
	DEPTH=DEPTH+1.0
	GOTO 5
C
	END
	
	
	
C*************************************************************************
	SUBROUTINE PHOZT(TIME,HMFELD,PM,DM,LFM,PT,DT,LFT,DEPTH)
C*************************************************************************
C
C	This subroutine calculates the depth of the photic zone inside of
C	the thermocline layer. Intervals of calculation are 1 meter.
C	The base of the photic zone is considered at that depth, where
C	1% of the surface PAR is attained. LPZ is the relative light intensity
C	at the base of the photic zone. Time enters in days, HMFELD in meters,
C	PM,DM,PT and DT are in mg Chla m-3, LFM and LFT are in Lith m-3, and
C	Depth is returned to the calling unit in meters.
C
	DOUBLE PRECISION TIME,HMFELD,PM,DM,LFM,PT,DT,LFT,DEPTH
	DOUBLE PRECISION Z,MIXL,HM,KEXR,KEXG,KEXB
	DOUBLE PRECISION VAR1,VAR2,VAR3,VAR4,VAR5,VAR6,VAR7
	DOUBLE PRECISION LPZ
C
	Z=0.0D0
C
	MIXL=HM(TIME,HMFELD)
C
	VAR1=EXP(-MIXL*KEXR(PM,DM))
	VAR2=EXP(-MIXL*KEXG(PM,DM,LFM))
	VAR3=EXP(-MIXL*KEXB(PM,DM,LFM))
	VAR4=(VAR1+VAR2+VAR3)/9.0
C
 5	VAR5=EXP(-Z*KEXR(PT,DT))
 	VAR6=EXP(-Z*KEXG(PT,DT,LFT))
	VAR7=EXP(-Z*KEXB(PT,DT,LFT))
C
	LPZ=VAR4*(VAR5+VAR6+VAR7)
C
C	Aphotic zone if the light intensity less
C	than 1 % of surface light intensity:
C
	IF (LPZ.LE.0.01) THEN
	  DEPTH=MIXL+Z
	  RETURN
	END IF
C
	Z=Z+1.0
	GOTO 5
C
	END
	
	
	



	
	
	
C***********************************************************************************************
	DOUBLE PRECISION FUNCTION MCM(TIME)
C***********************************************************************************************
C
C	Calculates the intrinsic grazing rates for coccolithophorids (MCM) in the mixed
C	layer according to the model of Taylor et al., 1990. MCM is returned to the
C	calling unit in day-1.
C
	DOUBLE PRECISION TIME,MCM
	DOUBLE PRECISION MCMMIN,MCMMAX
	DOUBLE PRECISION ARG,PI
C
	INCLUDE COCDAT.inc
C
	ARG=2.0*PI*(TIME-135.75)/365.0
C
	MCM=MCMMIN+0.5*(MCMMAX-MCMMIN)*(1.0+SIN(ARG))
C
	RETURN
	END
	
	
	
	
	
	
C***********************************************************************************************
	DOUBLE PRECISION FUNCTION MDM(TIME)
C***********************************************************************************************
C
C	Calculates the intrinsic grazing rates for diatoms (MDM) in the mixed
C	layer according to the model of Taylor et al., 1990. MDM is returned to the
C	calling unit in day-1.
C
	DOUBLE PRECISION TIME,MDM
	DOUBLE PRECISION MDMMIN,MDMMAX
	DOUBLE PRECISION ARG,PI
C
	INCLUDE COCDAT.inc
C
	ARG=2.0*PI*(TIME-135.75)/365.0
C
	MDM=MDMMIN+0.5*(MDMMAX-MDMMIN)*(1.0+SIN(ARG))
C
	RETURN
	END
	
	
	
C***********************************************************************************************
	DOUBLE PRECISION FUNCTION MCT(TIME)
C***********************************************************************************************
C
C	Calculates the intrinsic grazing rates for coccolithophorids (MCT) in the thermocline
C	layer according to the model of Taylor et al., 1990. MCT is returned to the
C	calling unit in day-1.
C
	DOUBLE PRECISION TIME,MCT
	DOUBLE PRECISION MCTMIN,MCTMAX
	DOUBLE PRECISION ARG,PI
C
	INCLUDE COCDAT.inc
C
	ARG=2.0*PI*(TIME-135.75)/365.0
C
	MCT=MCTMIN+0.5*(MCTMAX-MCTMIN)*(1.0+SIN(ARG))
C
	RETURN
	END




C***********************************************************************************************
	DOUBLE PRECISION FUNCTION MDT(TIME)
C***********************************************************************************************
C
C	Calculates the intrinsic grazing rates for diatoms (MDT) in the thermocline
C	layer according to the model of Taylor et al., 1990. MDT is returned to the
C	calling unit in day-1.
C
	DOUBLE PRECISION TIME,MDT
	DOUBLE PRECISION MDTMIN,MDTMAX
	DOUBLE PRECISION ARG,PI
C
	INCLUDE COCDAT.inc
C
	ARG=2.0*PI*(TIME-135.75)/365.0
C
	MDT=MDTMIN+0.5*(MDTMAX-MDTMIN)*(1.0+SIN(ARG))
C
	RETURN
	END





C***********************************************************************************************
	DOUBLE PRECISION FUNCTION ATMCOR(LAT)
C***********************************************************************************************
C
C	Calculates an atmospheric attenuation factor (reflection and absorption) 
C	for solar radiation at different latitudes (LAT). The input values are Ca 
C	(=solar radiation absorbed by clouds), Cr (=solar radiation reflected 
C	and scattered back to space by clouds), Aa (=solar radiation absorbed by 
C	dry air molecules, dust, and water vapour) and Ar (=solar radiation
C	reflected and scattered back to space by dry air molecules, dust and
C	water vapour). Ca, Cr, Aa and Ar were read from Sellers (1965) and are
C	stored in the external file ATM_COR. Output is the single value ATMCOR,
C	which indicates the fraction of solar radiation, which reaches the
C	Sea-surface (0.0˛ATRAN˛1.0). ATMCOR has the value of 1.0 at the top of
C	the atmosphere. LAT enters the function in degrees.
C
	INTEGER I
	DOUBLE PRECISION ATMCOR,LAT
	DOUBLE PRECISION ACOR(0:14)
	CHARACTER*15 HEADER
C
C	Reading the input data from file ATM_COR:
C
	OPEN(16,FILE='ATM_COR',STATUS='OLD')
	READ(16,*) HEADER
C
C	The vector ACOR(I) contains the latitude, and the values for Ca, Cr, Aa
C	and Ar at each latitude:
C
 250	READ(16,*,END=998) (ACOR(I),I=0,4)
 	IF (ACOR(0).EQ.LAT) THEN
	  CLOSE(16)
	  GOTO 300
	END IF
	GOTO 250
C
 300	ATMCOR=1.0-(ACOR(1)+ACOR(2)+ACOR(3)+ACOR(4))
C
	RETURN
C
 998	WRITE(9,*) 'End of ATM_COR reached'
 	STOP
C
	END
























































Include file COCDAT.inc (valid for COCDIA and PP programs)


C	Include file for COCDIA10.00.for and PP***.for
C	11.8.1994
C
	DOUBLE PRECISION DELTAT
	DOUBLE PRECISION MCMMIN,MCMMAX,MDMMIN,MDMMAX
	DOUBLE PRECISION MCTMIN,MCTMAX,MDTMIN,MDTMAX
	DOUBLE PRECISION VCOC,VDIA,K1S,K2S
	DOUBLE PRECISION GAMMAC,GAMMAD,EPSCOC,EPSDIA
	DOUBLE PRECISION TMAXM
	DOUBLE PRECISION TMINM,TB,NCHALF
	DOUBLE PRECISION PCMAX,AP,PDMAX,APD
	DOUBLE PRECISION KR0,KG0
	DOUBLE PRECISION LOSSIR,NDHALF,SHALF
	DOUBLE PRECISION HMIN,HMAX,HT,DLAT
	DOUBLE PRECISION CHLPC,CCPC,CCPCD,WLITH
C
C	The following variables are not ecological parameters
	DOUBLE PRECISION PI
C
C	DLAT: Latitude in degrees
	PARAMETER(DLAT=47.0)
	PARAMETER(DELTAT=1.0D0/24.0)
	PARAMETER(MCMMIN=0.04D0,MCMMAX=0.04D0)
	PARAMETER(MDMMIN=0.15D0,MDMMAX=0.15D0)
	PARAMETER(MCTMIN=0.04D0,MCTMAX=0.04D0)
	PARAMETER(MDTMIN=0.15D0,MDTMAX=0.15D0)
	PARAMETER(VCOC=2.7D-1)
C	K1S,K2S are summer (minimum) mixing coefficients, in m day-1
	PARAMETER(K1S=4.5D-1)
	PARAMETER(K2S=4.5D-1)
	PARAMETER(VDIA=1.5D0)
	PARAMETER(GAMMAC=5.03D-1)
	PARAMETER(EPSCOC=2.0D-1)
	PARAMETER(GAMMAD=1.424D0)
	PARAMETER(EPSDIA=2.0D-1)
	PARAMETER(TMAXM=17.0D0)
	PARAMETER(TMINM=9.0D0)
	PARAMETER(TB=8.0D0)
	PARAMETER(NCHALF=1.0D-1)
	PARAMETER(KR0=4.0D-1)
	PARAMETER(KG0=5.0D-2)
C	LOSSIR: Percent loss of infrared irradiation close to sea-surface
	PARAMETER(LOSSIR=0.5)
C	PCMAX: Maximum carbon assimilation rate of E.huxleyi at 15ˇC, in fmol C cell-1 hr-1
	PARAMETER(PCMAX=25.0D0)
C	AP: Specific absorption coefficient for chlorophyll in E.huxleyi, fmol C (cell hr)-1
	PARAMETER(AP=6.94D0)
C	PDMAX: Maximum carbon assimilation rate for diatoms, in fmol C cell-1 hr-1.
	PARAMETER(PDMAX=633.32D0)
C	APD: Specific absorption coefficient for chlorophyll in diatoms, calculated as
C	the best fit to approximate earlier Michaelis-Menten curve, in fmol C m-2 Ein-1 cell-1.
	PARAMETER(APD=800.0D0)
C	CCPCD: Organic carbon content per cell in diatoms, in fmol C cell-1.
	PARAMETER(CCPCD=7599.88D0)
	PARAMETER(NDHALF=3.0D-1)
	PARAMETER(SHALF=1.0D-1)
	PARAMETER(HMIN=20.0D0)
	PARAMETER(HT=40.0D0)
	PARAMETER(HMAX=3.5D2)
C	CHLPC: Inverse of chlorophyll content per cell in E. huxleyi, cell mgChla-1
	PARAMETER(CHLPC=6.11153D9)
C	CCPC: Organic carbon content of a single cell of E. huxleyi, in fmol C cell-1.
	PARAMETER(CCPC=544.89D0)
C	WLITH: Calcite mass of a single E. huxleyi coccolith, in grams
	PARAMETER(WLITH=3.46D-12)
	PARAMETER(PI=3.14159265D0)



































































	Alternative form of function GMAXE (TEMP)




        DOUBLE PRECISION FUNCTION GMAXE(TEMP)
C
C	Subroutine to calculate the maximum specific growth rate of
C	E. huxleyi as a function of temperature. Data were from Dave Lesley,
C	Roger Harris and Maureen Conte (Poster during fifth GEM Meeting at
C	Blagnac, September 1994). TEMP enters the subroutine in degrees
C	Celsius, and GMAXE is returned to the calling unit in day-1.
C	GMAXE stands for GMAX for E miliania huxleyi.
C
	INTEGER I,K,SI,EI
	DOUBLE PRECISION T(1:10),G(1:10)
	DOUBLE PRECISION GMAXE,TEMP,SLOPE
C
C	Reading temperatures of curve:
C
	DATA (T(I),I=1,10) /5.0D0,6.0D0,9.0D0,12.0D0,
     *  15.0D0,18.0D0,21.0D0,24.0D0,27.0D0,30.0D0/
C
C	Reading associated growth rates of curve:
C
	DATA (G(I),I=1,10) /0.0D0,0.144D0,0.243D0,
     *  0.466D0,0.592D0,0.773D0,1.024D0,0.981D0,0.935D0,0.0D0/
C
C	Calculating specific growth rates:
C	First, growth rates outside temperature tolerance range
C	for E. huxleyi are zero:
C
	IF ((TEMP.LT.T(1)).OR.(TEMP.GT.T(10))) THEN
	  GMAXE=0.0D0
	  RETURN
	END IF
C
C	Now, calculating growth rates if actual temperature is inside the
C	temperature tolerance range for E. huxleyi as function of TEMP:
C
	  DO 10, K=2,10
	    IF (TEMP.LE.T(K)) THEN
	      SI=K-1
	      EI=K
	      SLOPE=(G(EI)-G(SI))/(T(EI)-T(SI))
	      GMAXE=G(SI)+(TEMP-T(SI))*SLOPE
	      RETURN
	    END IF
 10	  CONTINUE
C
	END


















Input files to the COCDIA program




Input file INITCON.DAT with initial conditions

0.1,6.0D9,10.0,0.1,5.0,5.0
0.1,6.0D9,10.0,0.1,6.0,6.0
0.01,6.0D8,10.0,0.01,14.0,7.5
365.0,0.05,20









Input file MIXED_LAYER_DEPTHS

Explanation of columns:
There are 7 latitudinal stations.
The first column is the latitude, the next column is the starting mixed layer (day 0), the next 12 columns are monthly average mixed layer depths, the last column is the mixed layer depth at day 364.




35.O 111.5 134.4 180.3 135.0  97.5  63.2  35.6  26.9  30.0  36.1  46.5  60.0  90.0 111.5
40.0 126.7 160.9 229.4 215.4  90.0  45.0  30.0  30.0  30.0  38.9  45.0  64.5  93.0 126.7
45.0 119.0 198.9 358.8 418.8 277.8  90.0  28.4  30.0  24.4  42.0  58.3 112.5 159.1 119.0
47.0 212.0 283.0 425.0 500.0 420.0 122.0  29.0  22.5  25.4  42.7  67.5 135.0 190.9 212.0
50.0 367.5 418.8 521.4 500.0 412.5 129.0  34.5  28.7  31.0  45.0  78.9 150.0 350.0 367.5
55.0 425.0 464.3 542.9 500.0 441.7 314.3  75.0  30.0  39.9  45.0 105.0 200.0 366.7 425.0
60.0 614.0 655.6 738.9 778.6 650.0 400.0  55.0  30.0  30.0  53.4 130.0 200.0 481.2 614.0
 






























Input file SEA_SURFACE_ALBEDO 

In the file SEA_SURFACE_ALBEDO are the seasonal values for the albedo at a given latitude (i.e. the albedo at mid-time of each month). In the first column, the latitude is stored.




LAT 000 JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC 364
 90.0 0.00 0.00 0.00 0.00 0.20 0.15 0.12 0.10 0.10 0.18 0.00 0.00 0.00 0.00
 80.0 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08
 70.0 0.09 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.08 0.10 0.09
 60.0 0.11 0.10 0.10 0.08 0.08 0.07 0.07 0.07 0.07 0.08 0.08 0.10 0.11 0.11
 55.0 0.12 0.11 0.10 0.08 0.08 0.07 0.07 0.07 0.07 0.08 0.08 0.11 0.12 0.12
 50.0 0.12 0.11 0.10 0.08 0.07 0.06 0.06 0.06 0.07 0.07 0.08 0.11 0.12 0.12
 47.0 0.12 0.11 0.10 0.08 0.07 0.06 0.06 0.06 0.07 0.07 0.08 0.11 0.12 0.12
 45.0 0.12 0.11 0.10 0.08 0.07 0.06 0.06 0.06 0.07 0.07 0.08 0.11 0.12 0.12
 40.0 0.11 0.10 0.09 0.07 0.07 0.06 0.06 0.06 0.06 0.07 0.08 0.10 0.11 0.11
 35.0 0.10 0.10 0.08 0.07 0.07 0.06 0.06 0.06 0.06 0.07 0.08 0.09 0.10 0.10
 30.0 0.09 0.09 0.07 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07 0.08 0.09 0.09
 20.0 0.07 0.07 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07
 10.0 0.07 0.07 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07
  0.0 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06
-10.0 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06
-20.0 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.06 0.06 0.06 0.06 0.06 0.06
-30.0 0.06 0.06 0.06 0.07 0.08 0.09 0.09 0.09 0.08 0.07 0.07 0.06 0.06 0.06
-40.0 0.06 0.06 0.07 0.08 0.09 0.10 0.10 0.10 0.09 0.08 0.07 0.06 0.06 0.06
-50.0 0.06 0.06 0.07 0.08 0.09 0.11 0.11 0.11 0.09 0.08 0.07 0.06 0.06 0.06
-60.0 0.07 0.07 0.08 0.09 0.10 0.11 0.11 0.11 0.10 0.09 0.08 0.07 0.07 0.07
-70.0 0.55 0.55 0.55 0.55 0.80 0.80 0.80 0.80 0.80 0.80 0.55 0.55 0.55 0.55
-80.0 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80
-90.0 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80 0.80






































Input file ATM_COR

This file is called in function ATMCOR(LAT) fort he calculation of an atmospheric attenuation factor (reflection and absorption) for solar radiation at different latitudes (LAT).

Ca=solar radiation absorbed by clouds, Cr=solar radiation reflected and scattered back to space by clouds, Aa=solar radiation absorbed by dry air molecules, dust, and water vapour, and Ar=solar radiation reflected and scattered back to space by dry air molecules, dust and water vapour. Ca, Cr, Aa and Ar were read from Sellers (1965). 




Lat Ca Cr Aa Ar
90.0 0.030 0.273 0.121 0.091
80.0 0.029 0.297 0.116 0.087
70.0 0.033 0.303 0.118 0.092
60.0 0.033 0.286 0.115 0.088
55.0 0.030 0.273 0.119 0.083
50.0 0.027 0.259 0.123 0.077
47.0 0.026 0.250 0.125 0.074
45.0 0.026 0.245 0.127 0.072
40.0 0.024 0.230 0.131 0.067
35.0 0.021 0.215 0.138 0.062
30.0 0.018 0.199 0.145 0.057
20.0 0.020 0.198 0.158 0.050
10.0 0.022 0.221 0.164 0.044
0.0 0.025 0.236 0.168 0.040
-10.0 0.022 0.222 0.168 0.038
-20.0 0.020 0.211 0.158 0.050
-30.0 0.021 0.227 0.145 0.057
-40.0 0.028 0.276 0.126 0.067
-50.0 0.036 0.323 0.109 0.082
-60.0 0.033 0.348 0.105 0.088
-70.0 0.039 0.312 0.104 0.078
-80.0 0.029 0.217 0.087 0.029
-90.0 0.023 0.167 0.106 0.015





























Auxiliary programs:

Auxiliary program PP version 4.00 (PP400) for calculation of carbonate and Corg production




	PROGRAM PP
C
C	Version 4.00 created from PP300 on October 7, 1994.
C
C	Program to calculate the annual and vertically integrated carbonate and
C	Corg production by coccolithophorids and phytoplankton in the photic zone
C	simulated with program COCDIA (version 10.04).
C
C	Input data are required from files PLANKTON.DAT and PHYSICS.DAT, which are
C	both generated by COCDIA, and some parameters are read from the include file
C	COCDAT.inc. Output data are written to the external file PRODUCTION.DAT.
C
C	If COCDIA version 10.04 (daily light-dark cycle included) is used, then
C	it is necessary, that a record for every hour must be read from the input
C	files.
C
	INTEGER N
	DOUBLE PRECISION TIME,TIMEN,TIM,TIMN,PM,PMN,PT,PTN
	DOUBLE PRECISION DM,DMN,DT,DTN
	DOUBLE PRECISION HM,HMN,HP,HPN,MALPC,MALPCN,TALPC,TALPCN
	DOUBLE PRECISION MALPD,MALPDN,TALPD,TALPDN,PHICM,PHICMN
	DOUBLE PRECISION PHICT,PHICTN,PHIDM,PHIDMN,PHIDT,PHIDTN
	DOUBLE PRECISION GMXME,GMXMD,GMXTE,GMXTD
	DOUBLE PRECISION GMXMEN,GMXMDN,GMXTEN,GMXTDN
	DOUBLE PRECISION CALCM,CALCMN,CALCT,CALCTN
	DOUBLE PRECISION CALCIT,INCARB,ORGC,INCORG
C
	CHARACTER*36 H1PLANK
	CHARACTER*57 H2PLANK
	CHARACTER*96 H1PHYS
C
C	Open files "PLANKTON.DAT" (input), "PHYSICS.DAT" (input) and
C	"PRODUCTION.DAT" (output) and reading headers and initial pairs of
C	records from the input files:
C
	OPEN(15,FILE='PLANKTON.DAT',STATUS='OLD')
	OPEN(16,FILE='PHYSICS.DAT',STATUS='OLD')
	OPEN(17,FILE='PRODUCTION.DAT',STATUS='NEW')
C
C	Initialize counter N (for ouput to screen only):
C
	N=0
C
	READ(15,*) H1PLANK
	READ(15,*) H2PLANK
	READ(15,*) TIME,PM,PT,DM,DT
	READ(15,*) TIMEN,PMN,PTN,DMN,DTN
C
	READ(16,*) H1PLANK
	READ(16,*) H1PHYS
	READ(16,*) TIM,HM,HP,MALPC,TALPC,MALPD,TALPD,
     1  PHICM,PHICT,PHIDM,PHIDT,GMXME,GMXMD,GMXTE,
     2  GMXTD,CALCM,CALCT
     	READ(16,*) TIMN,HMN,HPN,MALPCN,TALPCN,MALPDN,TALPDN,
     3  PHICMN,PHICTN,PHIDMN,PHIDTN,GMXMEN,GMXMDN,GMXTEN,
     4  GMXTDN,CALCMN,CALCTN
C
C	Compute annual and vertically integrated carbonate and
C	organic carbon production:
C	Firstly, calculate initial increment (=CALCIT and ORGC):
C
	CALCIT=0.0D0
	ORGC=0.0D0
C
	CALL CARB(TIME,TIMEN,PM,PMN,PT,PTN,HM,HMN,HP,
     1            HPN,CALCM,CALCMN,CALCT,CALCTN,INCARB)
C
	CALCIT=CALCIT+INCARB
C
	CALL ORGCAR(TIME,TIMEN,PM,PMN,PT,PTN,DM,DMN,
     1  DT,DTN,HM,HMN,HP,HPN,MALPC,MALPCN,TALPC,TALPCN,
     2  MALPD,MALPDN,TALPD,TALPDN,PHICM,PHICMN,PHICT,PHICTN,
     3  PHIDM,PHIDMN,PHIDT,PHIDTN,GMXME,GMXMD,GMXMEN,
     4  GMXMDN,GMXTE,GMXTD,GMXTEN,GMXTDN,INCORG)
C
	ORGC=ORGC+INCORG
C
C	Writing initial increments to output file PRODUCTION.DAT:
C
	WRITE(17,*) 'Interval: TIME,TIMEN,CALCIT,ORGC'
	WRITE(17,*) TIME,TIMEN,CALCIT,ORGC
C
C	Writing message to screen:
C
	WRITE(9,*) '. . .Calculating. . .'
	WRITE(9,*) 'Time Interval = :',TIME,' - ',TIMEN
C
C	Now, applying "Strickleiter System": Replace all previous values (TIME,PM,DM,...)
C	by actual values (TIMEN,PMN,DMN...) and read next record from PLANKTON.DAT and
C	PHYSICS.DAT as actual values (TIMEN,PMN,DMN...):
C
 5	TIME=TIMEN
 	PM=PMN
	PT=PTN
	DM=DMN
	DT=DTN
	LAM=LAMN
	LAT=LATN
	LFM=LFMN
	LFT=LFTN
	NM=NMN
	NT=NTN
	SM=SMN
	ST=STN
	NB=NBN
	SB=SBN
	HM=HMN
	HP=HPN
	MALPC=MALPCN
	TALPC=TALPCN
	MALPD=MALPDN
	TALPD=TALPDN
	PHICM=PHICMN
	PHICT=PHICTN
	PHIDM=PHIDMN
	PHIDT=PHIDTN
	GMXME=GMXMEN
	GMXMD=GMXMDN
	GMXTE=GMXTEN
	GMXTD=GMXTDN
	CALCM=CALCMN
	CALCT=CALCTN
C
	READ(15,*,END=200) TIMEN,PMN,PTN,DMN,DTN
	READ(16,*,END=200) TIMN,HMN,HPN,MALPCN,TALPCN,
     1  MALPDN,TALPDN,PHICMN,PHICTN,PHIDMN,PHIDTN,
     2  GMXMEN,GMXMDN,GMXTEN,GMXTDN,CALCMN,CALCTN
C
C	Message to screen every 24 steps:
C
	IF (MOD(N,24).EQ.0.0) THEN
	  WRITE(9,*) 'Time Interval = :',TIME,' - ',TIMEN
	END IF
C
C	Calculate next increment in carbonate and organic carbon
C	production and add up:
C
	CALL CARB(TIME,TIMEN,PM,PMN,PT,PTN,HM,HMN,HP,
     1            HPN,CALCM,CALCMN,CALCT,CALCTN,INCARB)
C
	CALCIT=CALCIT+INCARB
C
	CALL ORGCAR(TIME,TIMEN,PM,PMN,PT,PTN,DM,DMN,
     1  DT,DTN,HM,HMN,HP,HPN,MALPC,MALPCN,TALPC,TALPCN,
     2  MALPD,MALPDN,TALPD,TALPDN,PHICM,PHICMN,PHICT,PHICTN,
     3  PHIDM,PHIDMN,PHIDT,PHIDTN,GMXME,GMXMD,GMXMEN,GMXMDN,
     4  GMXTE,GMXTD,GMXTEN,GMXTDN,INCORG)
C
	ORGC=ORGC+INCORG
C
C	Output to file:
C
	WRITE(17,*) TIME,TIMEN,CALCIT,ORGC
C
	N=N+1
C
	GOTO 5
C
 200	STOP
	END
	
	
	
	
	
C***********************************************************************************************
	SUBROUTINE CARB(TIME,TIMEN,PM,PMN,PT,PTN,HM,HMN,HP,
     1            HPN,CALCM,CALCMN,CALCT,CALCTN,INCARB)
C***********************************************************************************************
C
C	This subroutine calculates the increment in carbonate production by coccolithophorids
C	across the photic zone after the time increment DEL. Input values are PM,PMN,PT
C	PTN (mg Chla m-3), HM,HMN,HP,HPN (m), and CALCM,CALCMN,CALCT,CALCTN (Lith cell-1 day-1).
C	The increment in carbonate production is stored in INCARB and is returned to the calling
C	unit in gram calcite m-2.
C
	DOUBLE PRECISION TIME,TIMEN,PM,PMN,PT,PTN,HM,HMN
	DOUBLE PRECISION HP,HPN,CALCM,CALCMN,CALCT,CALCTN
	DOUBLE PRECISION INCARB,DEL,CHLPC,WLITH
	DOUBLE PRECISION XM,XMN,XT,XTN
C
	INCLUDE COCDAT.inc
C
C	CHLPC is the inverse of the chlorophyll content per cell in
C	E. huxleyi (in Cell mgChla-1).
C	WLITH: Calcite mass of a single E. huxleyi coccolith, in grams.
C
C	DEL is the time-interval between two subsequent data sets in the
C	input files:
C
	DEL=TIMEN-TIME
C
C	Convert mg Chla m-3 (PM,PMN,PT,PTN) into cell m-3 (XM,XMN,XT,XTN):
C
	XM=PM*CHLPC
	XMN=PMN*CHLPC
	XT=PT*CHLPC
	XTN=PTN*CHLPC
C
C	Carbonate production if the photic zone is entirely inside the mixed layer:
C
	IF (HPN.LE.HMN) THEN
	  INCARB=WLITH*DEL*(XM*CALCM*HP+XMN*CALCMN*HPN)/2.0
C
	ELSE
C	  Carbonate production if the photic zone reaches into the thermocline layer:
	  INCARB=WLITH*DEL*(XT*CALCT*(HP-HM)+
     1    XTN*CALCTN*(HPN-HMN)+
     2    XM*CALCM*HM+XMN*CALCMN*HMN)/2.0
	END IF
C
	RETURN
	END
	
	
	






	
	
C***********************************************************************************************
	SUBROUTINE ORGCAR(TIME,TIMEN,PM,PMN,PT,PTN,DM,DMN,
     1  DT,DTN,HM,HMN,HP,HPN,MALPC,MALPCN,TALPC,TALPCN,
     2  MALPD,MALPDN,TALPD,TALPDN,PHICM,PHICMN,PHICT,PHICTN,
     3  PHIDM,PHIDMN,PHIDT,PHIDTN,GMXME,GMXMD,GMXMEN,GMXMDN,
     4  GMXTE,GMXTD,GMXTEN,GMXTDN,INCORG)
C***********************************************************************************************
C
C	This subroutine calculates the increment in organic carbon production
C	by coccolithophorids and diatoms across the photic zone after the time
C	interval DEL. Phytoplankton concentrations (PM,PMN,DM,...) enter the
C	subroutine in mg Chla m-3, the specific growth rates are calculated in
C	day-1, INCORG is returned to the calling unit in gC m-2. The carbon to
C	chlorophyll ratio is 40 (already included in the factor of 0.02).
C
	DOUBLE PRECISION TIME,TIMEN,DEL,PM,PMN
	DOUBLE PRECISION PT,PTN,DM,DMN,DT,DTN
	DOUBLE PRECISION HM,HMN,HP,HPN,MALPC,MALPCN
	DOUBLE PRECISION TALPC,TALPCN,MALPD,MALPDN
	DOUBLE PRECISION TALPD,TALPDN,PHICM,PHICMN
	DOUBLE PRECISION PHICT,PHICTN,PHIDM,PHIDMN
	DOUBLE PRECISION PHIDT,PHIDTN,GMXME,GMXMD
	DOUBLE PRECISION GMXMEN,GMXMDN,GMXTE,GMXTD
	DOUBLE PRECISION GMXTEN,GMXTDN,INCORG
C
C	DEL is the time interval between two subsequent data sets
C	in the input files.
C
	DEL=TIMEN-TIME
C
C	Organic carbon production if the photic zone is
C	entirely inside the mixed layer:
C
	IF (HPN.LE.HMN) THEN
	  INCORG=0.02*DEL*(HP*(GMXME*PM*MALPC*PHICM+
     1	  GMXMD*DM*MALPD*PHIDM)+
     2    HPN*(GMXMEN*PMN*MALPCN*PHICMN+
     3    GMXMDN*DMN*MALPDN*PHIDMN))
	ELSE
	  INCORG=0.02*DEL*(HM*(GMXME*PM*MALPC*PHICM+
     4    GMXMD*DM*MALPD*PHIDM-
     5    GMXTE*PT*TALPC*PHICT-GMXTD*DT*TALPD*PHIDT)+
     6    HP*(GMXTE*PT*TALPC*PHICT+
     7    GMXTD*DT*TALPD*PHIDT)+
     8    HMN*(GMXMEN*PMN*MALPCN*PHICMN+
     9    GMXMDN*DMN*MALPDN*PHIDMN-
     1    GMXTEN*PTN*TALPCN*PHICTN-
     2    GMXTDN*DTN*TALPDN*PHIDTN)+
     3    HPN*(GMXTEN*PTN*TALPCN*PHICTN+
     4    GMXTDN*DTN*TALPDN*PHIDTN))
	END IF
C
	RETURN
	END
	




















Auxiliar program REDUCE (for data reduction)




	PROGRAM REDUCE
C
C	Version 2.00, 7. October, 1994
C
C	Program to reduce the hourly data file PRODUCTION.DAT
C	into a smaller data file PRODUCTION.RED.
C	Adapted to COCDIA Version 10.04 (different growth rate 
C	calculations for diatoms and coccolithophores, i.e. adjusted
C	to data file "PHYSICS.DAT")
C
	INTEGER N,NOUT
	REAL TIME,TIMEN,CALCIT,ORGC
	REAL PM,PT,DM,DT,LAM,LAT,LFM,LFT,NM,NT,SM,ST,NB,SB
C	
	REAL HM,HP,MALPHAC,TALPHAC,MALPHAD,TALPHAD
	REAL PHICM,PHICT,PHIDM,PHIDT,GMXME,GMXMD
	REAL GMXTE,GMXTD,CALCM,CALCT
C
	REAL ISURF,LIGHTM,LIGHTT,TEMPM,TEMPT,KEXRM,KEXGM,KEXBM
	REAL KEXRT,KEXGT,KEXBT
C
	CHARACTER*10 HEADER
C
	WRITE(9,*) '. . .Program to reduce hourly data sets. . .'
	WRITE(9,*) ' '
	WRITE(9,*) '. . .Degree of reduction:. . .'
	WRITE(9,*) '     NOUT= 1: No reduction'
	WRITE(9,*) '     NOUT=12: Output every 12 hours'
	WRITE(9,*) '     NOUT=24: Output every day, e.t.c'
	WRITE(9,*) ' '
	WRITE(9,*) '. . .Now enter NOUT. . .'
	READ(9,*) NOUT
C
C
C	First process PRODUCRTION.DAT:
C
	WRITE(9,*) '. . .Processing PRODUCTION.DAT. . .'
	OPEN(15,FILE='PRODUCTION.DAT',STATUS='OLD')
	OPEN(16,FILE='PRODUCTION.RED',STATUS='NEW')
C
C	Initialize counter:
C
	N=0
C
	READ(15,*) HEADER
C
C	Write header to output file:
	WRITE(16,*) 'Interval: TIME,TIMEN,CALCIT,ORGC'

 5	READ(15,*,END=999) TIME,TIMEN,CALCIT,ORGC
	IF (MOD(N,NOUT).EQ.0.0) THEN
	  WRITE(16,*) TIME,TIMEN,CALCIT,ORGC
	END IF
C
	N=N+1
	GOTO 5
C
 999	CLOSE(15)
	CLOSE(16)
C
C
C
C	Now process PLANKTON.DAT:
C
	WRITE(9,*) '. . .Processing PLANKTON.DAT. . .'
	OPEN(15,FILE='PLANKTON.DAT',STATUS='OLD')
	OPEN(16,FILE='PLANKTON.RED',STATUS='NEW')
C
C	Initialize counter:
C
	N=0
C
	READ(15,*) HEADER
	READ(15,*) HEADER
C
C	Write header to output file:
	WRITE(16,*) 'Time (Days),PM,PT,DM,DT,
     *  LAM,LAT,LFM,LFT,NM,NT,SM,ST,NB,SB'
C
 10	READ(15,*,END=888) TIME,PM,PT,DM,DT,
     *  LAM,LAT,LFM,LFT,NM,NT,SM,ST,NB,SB
	IF (MOD(N,NOUT).EQ.0.0) THEN
	  WRITE(16,*) TIME,PM,PT,DM,DT,
     *    LAM,LAT,LFM,LFT,NM,NT,SM,ST,NB,SB
	END IF
C
	N=N+1
	GOTO 10
C
 888	CLOSE(15)
	CLOSE(16)
C
C
C
C	Now process PHYSICS.DAT:
C
	WRITE(9,*) '. . .Processing PHYSICS.DAT. . .'
	OPEN(15,FILE='PHYSICS.DAT',STATUS='OLD')
	OPEN(16,FILE='PHYSICS.RED',STATUS='NEW')
C
C	Initialize counter:
C
	N=0
C
	READ(15,*) HEADER
	READ(15,*) HEADER
C
C	Write header to output file:
	WRITE(16,*) 'Time(days),HM,HP,MALPHAC,TALPHAC,
     *  MALPHAD,TALPHAD,PHICM,PHICT,PHIDM,PHIDT,
     *  GMXME,GMXMD,GMXTE,GMXTD,CALCM,CALCT'
C
 20	READ(15,*,END=777) TIME,HM,HP,MALPHAC,TALPHAC,
     *  MALPHAD,TALPHAD,PHICM,PHICT,PHIDM,PHIDT,
     *  GMXME,GMXMD,GMXTE,GMXTD,CALCM,CALCT
C
	IF (MOD(N,NOUT).EQ.0.0) THEN
	  WRITE(16,*) TIME,HM,HP,MALPHAC,TALPHAC,
     *    MALPHAD,TALPHAD,PHICM,PHICT,PHIDM,PHIDT,
     *    GMXME,GMXMD,GMXTE,GMXTD,CALCM,CALCT
	END IF
C
	N=N+1
	GOTO 20
C
 777	CLOSE(15)
	CLOSE(16)
C
C
C
C	Now process LIGHT.DAT:
C
	WRITE(9,*) '. . .Processing LIGHT.DAT. . .'
	OPEN(15,FILE='LIGHT.DAT',STATUS='OLD')
	OPEN(16,FILE='LIGHT.RED',STATUS='NEW')
C
C	Initialize counter:
C
	N=0
C
	READ(15,*) HEADER
C
C	Write header to output file:
	WRITE(16,*) 'TIME(days),ISURF,LIGHTM,LIGHTT,TEMPM,TEMPT,
     *  KEXRM,KEXGM,KEXBM,KEXRT,KEXGT,KEXBT'
C
 30	READ(15,*,END=666) TIME,ISURF,LIGHTM,LIGHTT,TEMPM,
     *  TEMPT,KEXRM,KEXGM,KEXBM,KEXRT,KEXGT,KEXBT
C
	IF (MOD(N,NOUT).EQ.0.0) THEN
	  WRITE(16,*) TIME,ISURF,LIGHTM,LIGHTT,TEMPM,
     *    TEMPT,KEXRM,KEXGM,KEXBM,KEXRT,KEXGT,KEXBT
	END IF
C
	N=N+1
	GOTO 30
C
 666	CLOSE(15)
	CLOSE(16)
C
	STOP
	END
	
	






































Auxiliary program TEST (for investigation of subroutines and functions)



	PROGRAM TEST
C
C	For testing subroutines and functions
C
	INTEGER N,NMAX,NOUT
	DOUBLE PRECISION DELTAT,LENGTH,TIME
	DOUBLE PRECISION TEMP,GMAXE,GROWTH
C
	INCLUDE COCDAT.inc
C
C	Determine length of simulation (DELTAT=1/24):
C
	WRITE(9,*) 'Enter length of simulation'
	WRITE(9,*) 'and output intervals (LENGTH,NOUT)'
	WRITE(9,*) ' '
	WRITE(9,*) 'Take care:'
	WRITE(9,*) 'NOUT= 1: output every hour'
	WRITE(9,*) 'NOUT=12: output every 12 hours'
	WRITE(9,*) 'NOUT=24: output every 24 hours'
	WRITE(9,*) '(i.e. every day at midnight)'
C
	READ(9,*) LENGTH,NOUT
C
C	Calculate maximum number of iterations (=NMAX)
C
	NMAX=INT(LENGTH/DELTAT)
C	Open output file 'TEST.DAT'
	OPEN(15,FILE='TEST.DAT',STATUS='NEW')
	WRITE(15,*) 'TIME(days),TEMP,GMAXE'
C
	WRITE(9,*) '. . .Calculating. . .'
C
C	Initialization of iteration:
	N=0
C
C	Iteration:
 5	TIME=N*DELTAT
C****
	TEMP=TIME
	GROWTH=GMAXE(TEMP)
C****
C
C	Output:
C
	IF (MOD(N,NOUT).EQ.0.0) THEN
	  WRITE(15,*) TIME,TEMP,GROWTH
	END IF
C
C	Stop iteration:
C
	IF (N+1.GE.NMAX) THEN
	  WRITE(9,*) 'One file "TEST.DAT" written'
	  WRITE(9,*) 'Enter any key to continue'
	  PAUSE 999
 999	  CONTINUE
	  STOP
	END IF
C
C	Calculate the next time-step:
C
	N=N+1
	GOTO 5
C
	END
	DOUBLE PRECISION FUNCTION GMAXE(TEMP)
C
C	Soubroutine to calculate the maximum specific growth rate of
C	E. huxleyi as a function of temperature. Data were from Dave Lesley,
C	Roger Harris and Maureen Conte (Poster during fifth GEM Meeting at
C	Blagnac, September 1994). TEMP enters the subroutine in degrees
C	Celsius, and GMAX is returned to the calling unit in day-1.
C	GMAXE stands for GMAX for E miliania huxleyi.
C
	INTEGER I,K,SI,EI
	DOUBLE PRECISION T(1:10),G(1:10)
	DOUBLE PRECISION GMAXE,TEMP,SLOPE
C
C	Reading temperatures of curve:
C
	DATA (T(I),I=1,10) /5.0D0,6.0D0,9.0D0,12.0D0,
     *  15.0D0,18.0D0,21.0D0,24.0D0,27.0D0,30.0D0/
C
C	Reading associated growth rates of curve:
C
	DATA (G(I),I=1,10) /0.0D0,0.144D0,0.243D0,
     *  0.466D0,0.592D0,0.773D0,1.024D0,0.981D0,0.935D0,0.0D0/
C
C	Calculating specific growth rates:
C	First, growth rates outside temperature tolerance range
C	for E. huxleyi are zero:
C
	IF ((TEMP.LT.T(1)).OR.(TEMP.GT.T(10))) THEN
	  GMAXE=0.0D0
	  RETURN
	END IF
C
C	Now, calculating growth rates if actual temperature is inside the
C	temperature tolerance range for E. huxleyi as function of TEMP:
C
	  DO 10, K=2,10
	    IF (TEMP.LE.T(K)) THEN
	      SI=K-1
	      EI=K
	      SLOPE=(G(EI)-G(SI))/(T(EI)-T(SI))
	      GMAXE=G(SI)+(TEMP-T(SI))*SLOPE
	      RETURN
	    END IF
 10	  CONTINUE
C
	END




1