-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFlexEFT.f90
executable file
·1650 lines (1398 loc) · 53.1 KB
/
FlexEFT.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
PROGRAM main
implicit none
character(LEN=2 ), parameter :: Stn = 'K2'
! Input files
character(LEN=20), parameter :: Tempfile = trim(Stn)//'_Temp.dat', &
PARfile = trim(Stn)//'_par.dat' , &
Aksfile = trim(Stn)//'_Aks.dat' , &
NO3file = trim(Stn)//'_NO3.dat'
! !INPUT PARAMETERS for diffusion:
integer,parameter :: Dirichlet = 0
integer,parameter :: Neumann = 1
integer, parameter :: nlev = 40, & ! Number of vertical layers
NVAR = 7, & ! Number of total state variables
! Indices for state variables
iNO3 = 1, iPHY = 2, iZOO = 3, iDET = 4, &
iPMU = 5, iVAR = 6, iCHL = 7, &
NVsinkterms = 2, &
Windex(NVsinkterms) = (/iPHY, iDET/), &
! Number of vertical points of NO3 profile
N_NO3 = 14, &
N_Temp = 24, &
N_Aks = 40, &
Nout = 28, & ! Number of output variables
oTemp = 1, &
oPAR = 2, &
oAks = 3, &
oNO3 = 1, &
oPHY = 2, &
oZOO = 3, &
oDET = 4, &
oPMU = 5, &
oVAR = 6, &
oCHL = 7, &
oFER = 8, &
oTheta = 9, &
oQN = 10, &
omuNet = 11, &
oGraz = 12, &
oZ2N = 13, &
oD2N = 14, &
odmudl = 15, &
od2mu = 16, &
od2gdl = 17, &
ow_p = 18, &
! Diffusion fluxes
oD_NO3 = 19, &
oD_PHY = 20, &
oD_ZOO = 21, &
oD_DET = 22, &
oD_PMU = 23, &
oD_VAR = 24, &
oD_CHL = 25
! Output indices:
character(LEN=5), parameter :: Labelout(Nout) &
= (/'Temp ','PAR ','Aks ','NO3 ','PHY ','ZOO ','DET ','PMU ',&
'VAR ','CHL ','FER ','Theta','QN ','muNet','Graz ','Z2N ',&
'D2N ','dmudl','d2mu ','d2gdl','w_p ', &
'D_NO3','D_PHY','D_ZOO','D_DET','D_PMU','D_VAR','D_CHL'/)
real(8), parameter :: hmax = 250., & ! Total water depth
thetaS= 2. , & ! surface stretching parameter
dtsec = 60. , & ! time step in seconds
d_per_s = 86400.,& ! how many seconds in one day
dtdays= dtsec/d_per_s, & !the fraction of a time step in one day
wDET = 2d0 ! Sinking rate of Detritus (m/d)
integer, parameter :: y_per_s = INT(d_per_s*360), &!how many seconds of one year
nsave= INT(d_per_s)/INT(dtsec), &! Timesteps to save
maxD = 360*5, & ! Number of days to run
maxN = maxD*INT(d_per_s)/INT(dtsec) ! Number of time steps
real(8),parameter :: DefaultRelaxTau(nlev)= 1.d15, cnpar = 0.6
real(8) :: Z_r(1:nlev), Z_w(0:nlev), Hz(nlev), & ! Grid variables
Vars(NVar,nlev), & ! State variables
NewVars(NVar,nlev) , & ! New calculated state variables
obs_NO3(N_NO3 , 2), &
obs_PAR(1 ,13), &
obs_Temp(N_Temp,13), &
obs_Aks(N_Aks ,13), &
obs_time(12), & ! observation time indices
time, I_0, PAR(nlev), Temp(nlev), Aks(0:nlev), &
NO3(nlev,1), &
VTemp(nlev,12), VAks(0:nlev,12), Varout(100,nlev), &
! Sinking rate
ww(0:nlev,NVsinkterms), Lsour(nlev), Qsour(nlev)
integer :: i, it, k, j, time_day
character(LEN=20) :: outfile
! Setup grid
call setup_grid(thetaS, nlev, hmax, Z_r, Z_w, Hz)
! Initialize state variables
call Readcsv(NO3file, N_NO3, 2, obs_NO3)
! Interpolate initial NO3:
call gridinterpol(N_NO3,1,obs_NO3(:,1),obs_NO3(:,2), &
nlev, Z_r, NO3(:,1))
! Initialize initial NO3:
Vars(iNO3,:) = NO3(:,1)
! Initialize other variables:
do i = 2,NVAR
do k = 1,nlev
Vars(i, k) = 0.1
enddo
enddo
! Initialize output data:
Varout(:,:) = 0.0
do i = 1,NVAR
do k = 1,nlev
Varout(i,k)= Vars(i,k)
enddo
enddo
! Give obs. time indices in seconds:
do i = 1, 12
obs_time(i) = (float(i-1)+0.5)*30.*d_per_s
enddo
! Read external PAR data:
call Readcsv(PARfile, 1, 13,obs_PAR)
! Read external Temp data:
call Readcsv(Tempfile,N_Temp, 13,obs_Temp)
! Interpolate external Temp data:
call gridinterpol(N_Temp,12,obs_Temp(:,1),obs_Temp(:,2:13), &
nlev, Z_r, VTemp(:,:))
! Interpolate external Aks data:
call Readcsv(Aksfile, N_Aks , 13,obs_Aks)
! Interpolate external Aks data:
call gridinterpol(N_Aks,12,obs_Aks(:,1), obs_Aks(:,2:13), &
nlev+1, Z_w, VAks(:,:))
! Create output files:
do i = 1, Nout
outfile = trim(Labelout(i))//'.out'
! Create data out file
open (9+i, file = outfile, status = 'replace')
if (i .eq. oAks) then
write(9+i, 3) 'Timestep','Days',Z_w
else
write(9+i, 1) 'Timestep','Days',Z_r
endif
enddo
1 format(A10,1x,A7, <nlev >(2x,F12.5))
2 format(I10,1x,I7, <nlev >(2x,F12.5))
3 format(A10,1x,A7, <nlev+1>(2x,F12.5))
4 format(I10,1x,I7, <nlev+1>(2x,F12.5))
! For each time step, read in external environmental data
write(6,*) 'START TIME STEPPING'
DO it = 1, maxN+1
! Calculate current timing (zero is starting time):
time = float(it-1)*dtsec
! Interpolate Temp data throughout the water column:
call time_interp(int(time), 12, nlev, obs_time, VTemp, Temp)
! Interpolate temporal PAR data:
call time_interp(int(time), 12,1, obs_time,obs_PAR, I_0)
! Convert the unit of par to W m-2:
I_0 = I_0/0.4
! Calculate light field:
call Calculate_PAR(I_0, nlev, Hz, Vars(iCHL,:), PAR)
! Interpolate Aks data throughout the water column:
call time_interp(int(time), 12, nlev+1, obs_time, VAks, Aks)
! Save data to output files:
if (mod(it, nsave) .eq. 1) then
! Calculate model time in days:
time_day=int(time/d_per_s)
! Save Temp data into the output file:
write(9+oTemp, 2) it, time_day, Temp
! Save PAR data into the output file:
write(9+oPAR, 2) it, time_day, PAR
! Save Aks data into the output file:
write(9+oAks, 4) it, time_day, Aks
! Save state variables and diagnostics:
do i=4,Nout
write(9+i, 2) it, time_day, Varout(i-3,:)
enddo
endif !==> End of saving results
! Biological rhs:
do k = 1,nlev
call FlexEFT(Temp(k),PAR(k),dtdays,NVAR,Vars(:,k),Varout(:,k))
! Pass the new state variables to Vars
do j = 1,NVAR
Vars(j,k)=Varout(j,k)
enddo
enddo
! Diffusion:
Lsour(:) = 0.
Qsour(:) = 0.
do j = 1,NVAR
call diff_center(nlev,dtsec,cnpar,1,Hz,Neumann,Neumann,&
0.,0.,Aks,Lsour,Qsour,DefaultRelaxTau,Vars(j,:),Vars(j,:),NewVars(j,:))
! Save diffusion fluxes (normalized to per day)
Varout(18+j,:) = (NewVars(j,:) - Vars(j,:))/dtdays
! Update the state variables:
Vars(j,:)=NewVars(j,:)
enddo
! Sinking:
! Initialize sinking rate:
ww(0, :) = 0.0
ww(nlev,:) = 0.0
do k = 1,nlev-1
ww(k,1) = Varout(ow_p,k) !Phytoplankton sinking rate
ww(k,2) = -wDET !Detritus sinking rate
enddo
do j = 1,NVsinkterms
call adv_center(nlev,dtdays,Hz,Hz,ww(:,j), &
1,1,0.,0.,6,1,Vars(Windex(j),:))
enddo
! Check whether the values are valid:
do j = 1,NVAR
do k=1,nlev
if( (Vars(j,k) .ne. Vars(j,k)) .OR. (Vars(j,k) .le. 0.) ) then
write(6,*) 'At time step ',it
write(6,*) 'The variable ',trim(Labelout(j+3)), 'is invalid at depth ',Z_r(k)
stop
endif
enddo
enddo
ENDDO ! ==> End of time stepping
! Close files:
do i = 1, Nout
close (unit=9+i)
enddo
END PROGRAM main
!========================================================
! !ROUTINE: Interpolate from time-series observation to model time
subroutine time_interp(time, N, nrow, obs_time, obs_data, mod_data)
! time: model time
! N: number of time-series observations in observational data
! nrow: number of vertical points in observation data
! obs_time: time of observational data, should be a vector of length N
! obs_data: observational data, should be a vector of length N
! mod_data: model data as output
implicit none
integer, intent(in) :: N, time, nrow
real(8), intent(in) :: obs_time(N), obs_data(nrow,N)
real(8), intent(out) :: mod_data(nrow)
integer, parameter :: y_per_s = 31104000
integer :: i,j
real(8) :: timeMOD, rat
! Get time of the year (mod):
timeMOD = float(mod(time,y_per_s))
! Deal with the case between two years:
IF ((timeMOD .lt. obs_time(1))) THEN
rat = (timeMOD -obs_time(N)+float(y_per_s)) &
/ (obs_time(1)-obs_time(N)+float(y_per_s))
do i=1,nrow
mod_data(i)=obs_data(i,N)*(1d0-rat)+obs_data(i,1)*rat
enddo
ELSEIF (timeMOD .ge. obs_time(N)) THEN
rat = (timeMOD -obs_time(N)) &
/ (obs_time(1)-obs_time(N)+float(y_per_s))
do i=1,nrow
mod_data(i)=obs_data(i,N)*(1d0-rat)+obs_data(i,1)*rat
enddo
ELSE
do j = 1, N-1
if ((timeMOD.lt.obs_time(j+1)) .and. (timeMOD.GE.obs_time(j))) then
rat=(timeMOD-obs_time(j))/(obs_time(j+1)-obs_time(j))
do i=1,nrow
mod_data(i)=obs_data(i,j)*(1d0-rat)+obs_data(i,j+1)*rat
enddo
exit
endif
enddo
ENDIF
return
end subroutine
!========================================================
! !ROUTINE: Interpolate from observation space to model grid
!
! !INTERFACE:
subroutine gridinterpol(N,cols,obs_z,obs_prof,nlev,model_z,model_prof)
!
! !DESCRIPTION:
!
! This is a utility subroutine in which observational data, which might
! be given on an arbitrary, but structured grid, are linearly interpolated and
! extrapolated to the actual (moving) model grid.
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!N : number of vertical layers in observation data
!cols : number of profiles in obs. data
integer, intent(in) :: N,cols
REAL(8), intent(in) :: obs_z(N),obs_prof(N,cols)
integer, intent(in) :: nlev
REAL(8), intent(in) :: model_z(nlev)
!
! !OUTPUT PARAMETERS:
REAL(8), intent(out) :: model_prof(nlev,cols)
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
!
!EOP
!
! !LOCAL VARIABLES:
integer :: i,j,ii
REAL(8) :: rat
!-----------------------------------------------------------------------
!BOC
! Set surface values to uppermost input value
do i=nlev,1,-1
if(model_z(i) .ge. obs_z(N)) then
do j=1,cols
model_prof(i,j) = obs_prof(N,j)
end do
end if
end do
! Set bottom values to lowest input value
do i=1,nlev
if(model_z(i) .le. obs_z(1)) then
do j=1,cols
model_prof(i,j) = obs_prof(1,j)
end do
end if
end do
! Interpolate inner values linearly
do i=1,nlev
!write(*,*) 'i = ',i
if ((model_z(i) .lt. obs_z(N)) .and. (model_z(i) .gt. obs_z(1))) then
ii=0
224 ii=ii+1
if (obs_z(ii) .le. model_z(i)) goto 224
! Debug:
!write(*,*) 'ii = ', ii
!write(*,*) 'model_z = ',model_z(i)
!write(*,*) ' obs_z = ', obs_z(ii)
!write(*,*) ' obs_z(ii-1) = ', obs_z(ii-1)
rat=(model_z(i)-obs_z(ii-1))/(obs_z(ii)-obs_z(ii-1))
do j=1,cols
model_prof(i,j)=(1-rat)*obs_prof(ii-1,j)+rat*obs_prof(ii,j)
!write(*,*) 'obs_prof(ii-1) = ',obs_prof(ii-1,j)
!write(*,*) 'obs_prof(ii) = ',obs_prof(ii ,j)
!write(*,*) 'rat = ',rat
!write(*,*) 'model_prof(i) = ',model_prof(i,j)
end do
end if
end do
return
end subroutine
!EOC
!========================================================
subroutine setup_grid(thetaS, nlev, hmax, Z_r, Z_w, Hz)
implicit none
integer, intent(in) :: nlev
real(8), intent(in) :: thetaS, hmax
real(8), intent(out) :: Z_r(1:nlev), Z_w(0:nlev), Hz(1:nlev)
! Local scratch variable:
integer :: i
real(8) :: sc_r, C_sig
Z_w(0) = -hmax
!Following Song and Haidvogel (1994). sinh is the hyperbolic sin function
do i = 1,nlev
sc_r = (float(i-nlev) - 0.5)/float(nlev)
C_sig = sinh(thetaS*sc_r)/sinh(thetaS) ! -1 < C_sig < 0
Z_r(i) = C_sig*hmax
sc_r = (float(i-nlev))/float(nlev)
C_sig = sinh(thetaS*sc_r)/sinh(thetaS) ! -1 < C_sig < 0
Z_w(i) = C_sig*hmax
Hz(i) = Z_w(i) - Z_w(i-1)
enddo
return
end subroutine
!========================================================
SUBROUTINE FlexEFT(tC,par,dtdays,NVAR,Varin,Varout)
implicit none
!INPUT PARAMETERS:
real(8), intent(in) :: tC,par,dtdays
integer, intent(in) :: NVAR
real(8), intent(in) :: Varin(NVAR)
integer, parameter :: Nout = 100
real(8), intent(out) :: Varout(Nout)
!LOCAL VARIABLES of phytoplankton:
real(8) :: NO3,PHY,ZOO,DET,PMUPHY,VARPHY,NO31,PHY1,ZOO1,DET1, &
PMUPHY1,VARPHY1,CHL,CHL1 !State variables
real :: tf_p,tf_z,I_zero,larg !Environmental variables
real(8) :: PMU,VAR,PMU1,VAR1,X,B,dXdl,dBdl,KFe,Fe,dmu0hatdl, d2mu0hatdl2
real(8) :: V0hat,Kn,Lmin,A0N,A0hat,dA0hatdl,d2A0hatdl2,fA
real(8) :: VNhat,dVNhatdl,d2VNhatdl2, fN, d2fNdl2 ! Nutrient uptake variables
real(8) :: mu0hat,muIhat,mu0hatSI,dmu0hatSIdl,d2mu0hatSIdl2
real(8) :: dmuIhatdl,d2muIhatdl2! Growth rate variables
real(8) :: fV,dfVdl,d2fVdl2 !fV
real(8) :: ZINT,dZINdl,d2ZINdl2 !ZINT
real(8) :: aI,SI,dSIdl,d2SIdl2 !Light dependent component
real(8) :: RMchl,Theta,ThetaHat,dThetaHatdl,d2ThetaHatdl2 !Chl related variables
real(8) :: QN,Qs,dQNdl,d2QNdl2 ! cell quota related variables
real(8) :: larg1,w_p,dmu0hat_aIdl,d2mu0hat_aIdl2,dlargdl, &
d2largdl2,W_Y,dWYYdl,daI_mu0hatdl,d2aI_mu0hatdl2, &
d2wpdl2
real(8) :: dV0hatdl,d2V0hatdl2,muNet,dmuNetdl,d2muNetdl2,ThetaAvg,QNavg,dwdl, &
d2wdl2,w_pAvg,dgdlbar,d2gdl2bar
real(8),parameter :: thetamin = 0.01, DETmin = 1D-2,DETmax = 1D5 &
,QNmin = 0.01,fVmin = 0.01, Lref = 0.5, eps = 1D-4 &
,PMUmax= 20., VARmax = 100., Penfac = 1.0, pi = 3.14159265358979 &
,Ep = 0.4, alphaQ = -0.18,alphaI = -0.13, alphamu = 0. &
,betamu= 0. , Femin = 0.02, mu0 = 5.0, aI0 = 0.25 &
,A0N0 = 5.0, alphaA = -0.3, V0N = 5.0, alphaV = 0.2 &
,K0N = 0.17, alphaK = 0.27, K0Fe = 0.8,alphaFe = 0.14 &
,alphaG = 1.1, zetaChl= 0.8 , zetaN = 0.6,RMchl0 = 0.1 &
,mp = 0.0, Q0N = 0.1, w_p0 = 0.0, alphaW = 0.0, &
Wmax = 50.0, gmax = 1.0, kp = 0.5, &
rdn = 0.2, Ez = 0.6, mz = 0.05
integer, parameter :: nutrient_uptake = 2, grazing_formulation = 3
real(8), external :: temp, WAPR, ScaleTrait, PenAff, grazing
!Declarations of zooplankton:
real(8) :: Cf,aG,INGES,RES,EGES,gbar
real(8),parameter :: unass = 0.24, GGE = 0.3333
integer :: i,j,gform
integer, parameter :: iNO3=1,iPHY=2,iZOO=3,iDET=4,iPMU=5,iVAR=6,iCHL=7,&
iFER=8, &
oNO3 = 1, &
oPHY = 2, &
oZOO = 3, &
oDET = 4, &
oPMU = 5, &
oVAR = 6, &
oCHL = 7, &
oFER = 8, &
oTheta = 9, &
oQN = 10, &
omuNet = 11, &
oGraz = 12, &
oZ2N = 13, &
oD2N = 14, &
odmudl = 15, &
od2mu = 16, &
od2gdl = 17, &
ow_p = 18
real(8) :: pp(4,4), PP_pn,Zmort
logical, parameter :: do_IRON = .false.
!-----------------------------------------------------------------------
!BOC
! Retrieve current (local) state variable values.
NO3 = Varin(iNO3)
PHY = Varin(iPHY)
ZOO = Varin(iZOO)
DET = Varin(iDET)
CHL = Varin(iCHL)
PMUPHY = Varin(iPMU)
VARPHY = Varin(iVAR)
Varout(:) = 0.
!All rates have been multiplied by dtdays to get the real rate correponding to the actual time step
tf_p = temp(Ep,tC)
DET = max(DET,eps)
PHY = max(PHY,eps)
PMUPHY = max(PMUPHY,eps)
VARPHY = max(VARPHY,eps)
PMU = PMUPHY/PHY
VAR = VARPHY/PHY
! Fe related:
if (do_IRON .eq. .true.) then
Fe = Varin(iFER)
Fe = max(Fe,Femin) !Dissolved Fe concentration
KFe = ScaleTrait(PMU, K0Fe,alphaFe) !The half saturation constant of Fe at average size
endif
Qs = ScaleTrait(PMU, Q0N, alphaQ)/2.
mu0hat = dtdays*tf_p*mu0*exp(alphamu*PMU + betamu*PMU*PMU)
X = 0.
if (do_IRON .eq. .true.) then
mu0hat = mu0hat * Fe/(Fe + KFe)
X = alphaFe*KFe/(Fe + KFe)
endif
dmu0hatdl = mu0hat*(alphamu + 2D0*betamu*PMU - X)
d2mu0hatdl2= dmu0hatdl*(alphamu + 2D0*betamu*PMU-X) + mu0hat*2D0*betamu
if (do_IRON .eq. .true.) then
d2mu0hatdl2 = d2mu0hatdl2 - mu0hat*alphaFe*Fe*X/(Fe+KFe)
endif
! Iron limits nitrogen uptake
V0hat = ScaleTrait(PMU, dtdays*tf_p*V0N, alphaV)
if (do_IRON .eq. .true.) V0hat = V0hat*Fe/(Fe + KFe)
dV0hatdl = V0hat*(alphaV- X)
d2V0hatdl2 = dV0hatdl*(alphaV - X)
if (do_Iron .eq. .true.) then
d2V0hatdl2 = d2V0hatdl2 - V0hat*alphaFe*Fe*X/(Fe+KFe)
endif
! Initial slope of P-I curve
aI = ScaleTrait(PMU, dtdays*aI0, alphaI)
! Cost of photosynthesis
RMchl = tf_p*RMchl0*dtdays
! Threshold irradiance and RMchl is set temperature dependent
I_zero = zetaChl*RMchl/aI
!Define VNhat: the saturation function of ambient nutrient concentration
selectcase(nutrient_uptake)
! case 1: Classic Michaelis Menton
case(1)
! Potential maximal nutrient-limited uptake
! Half-saturation constant of nitrate uptake
Kn = ScaleTrait(PMU,K0N, alphaK)
VNhat = V0hat*NO3/(NO3 + Kn)
dVNhatdl = -VNhat*alphaK*Kn/(NO3+Kn) + NO3/(NO3 + Kn)*dV0hatdl
d2VNhatdl2 = -alphaK*(VNhat * alphaK * NO3/(NO3+Kn)**2*Kn &
+ Kn/(NO3 + Kn)*dVNhatdl)+ NO3/(NO3+ Kn)*d2V0hatdl2 &
- dV0hatdl*NO3/(NO3 + Kn)**2*alphaK*Kn
! case 2: optimal uptake based on Pahlow (2005) and Smith et al. (2009)
case(2)
Lmin = log(Lref**3/6.*pi) + log(1.+Penfac)/( Penfac*alphaA)
A0N = dtdays*tf_p*A0N0
A0hat = PenAff(PMU, alphaA, Penfac,Lmin) * ScaleTrait(PMU, A0N, alphaA)
A0hat = max(A0hat,eps)
dA0hatdl = alphaA*A0hat-A0N*exp(PMU*alphaA)*Penfac*alphaA &
* exp(Penfac*alphaA *(PMU-Lmin))
d2A0hatdl2 = alphaA*dA0hatdl &
- Penfac*alphaA*exp(alphaA*((1.+Penfac)*PMU-Penfac*Lmin)) &
* (dA0hatdl + A0hat*alphaA*(1.+Penfac))
!Define fA
fA = 1D0/( 1D0 + sqrt(A0hat * NO3/V0hat) )
VNhat = (1D0-fA)*V0hat*fA*A0hat*NO3/((1D0-fA)*V0hat + fA*A0hat*NO3)
!X: temporary variable
X = V0hat/A0hat + 2.*sqrt(V0hat*NO3/A0hat) + NO3
!B: d(V0/A0)dl
B = dV0hatdl/A0hat - V0hat/A0hat**2*dA0hatdl
dXdl = B*(1. + sqrt(NO3 * A0hat / V0hat))
dBdl = d2V0hatdl2/A0hat - dV0hatdl*dA0hatdl/A0hat**2 &
- (V0hat/A0hat**2*d2A0hatdl2 &
+ dA0hatdl*(dV0hatdl/A0hat**2 - 2.*V0hat*dA0hatdl/A0hat**3))
dVNhatdl = NO3*(dV0hatdl/X-V0hat/X**2*B*(1.+ sqrt(NO3*A0hat/V0hat) ) )
d2VNhatdl2 = NO3*(d2V0hatdl2/X - dV0hatdl*dXdl/X**2 &
- (V0hat/X**2*B*(-sqrt(NO3)*.5*(A0hat/V0hat) &
* sqrt(A0hat/V0hat) * B) &
+ B*(1. + sqrt(NO3*A0hat/V0hat) ) &
* (dV0hatdl/X**2 - 2.*V0hat*dXdl/X**3) &
+ V0hat/X**2*(1.+sqrt(NO3*A0hat/V0hat))*dBdl))
case default
write(*,*) 'Error: Incorrect option for nutrient uptake!'
ENDSELECT
! Calculate thetahat (optimal g Chl/mol C for the chloroplast under nutrient replete conditions)
! Only calculate within the euphotic zone, otherwise many numerical problems.
if( par .gt. I_zero ) then
larg1 = exp(1. + min(aI*par/(mu0hat*zetaChl),600.))
larg = (1. + RMchl/mu0hat)*larg1
dmu0hat_aIdl = (dmu0hatdl - alphaI*mu0hat)/aI
d2mu0hat_aIdl2 = d2mu0hatdl2/aI -alphaI/aI*dmu0hatdl-aI*dmu0hat_aIdl
daI_mu0hatdl = -(aI/mu0hat)**2*dmu0hat_aIdl
d2aI_mu0hatdl2 = -((aI/mu0hat)**2*d2mu0hat_aIdl2 &
- 2./(mu0hat/aI)**3*(dmu0hat_aIdl)**2)
dlargdl = -RMchl*larg1/(mu0hat**2)*dmu0hatdl &
+ (1.+RMchl/mu0hat)*larg1 * par/zetaChl*daI_mu0hatdl
d2largdl2 = -RMchl*(larg1*mu0hat**(-2)*d2mu0hatdl2 &
+ larg1*par/zetaChl*daI_mu0hatdl*mu0hat**(-2)*dmu0hatdl &
+ larg1*dmu0hatdl*(-2.*mu0hat**(-3)*dmu0hatdl)) &
+ par/zetaChl*((1+RMchl/mu0hat)*larg1*d2aI_mu0hatdl2 &
+ (1.+RMchl/mu0hat)*larg1*par/zetaChl*daI_mu0hatdl*daI_mu0hatdl &
+ RMchl*(-mu0hat**(-2)*dmu0hatdl)*larg1*daI_mu0hatdl)
W_Y = WAPR(larg,0,0)
ThetaHat = 1.0/zetaChl + (1.0- W_Y)*mu0hat/(aI * par)
ThetaHat = max(ThetaHat,thetamin)
dThetaHatdl = 1./par*(-W_Y/larg/(1.+W_Y)*dlargdl*mu0hat/aI &
+ (1.-W_Y)*dmu0hat_aIdl)
dWYYdl = dlargdl*(-W_Y**2/larg**2/(1.+W_Y)**3 &
- W_Y/larg**2/(1.+W_Y) + W_Y/(larg*(1.+W_Y))**2)
d2ThetaHatdl2 = 1/par*(-(W_Y/larg/(1.+W_Y)*dlargdl*dmu0hat_aIdl &
+ W_Y/larg/(1.+W_Y)*d2largdl2*mu0hat/aI &
+ dWYYdl*dlargdl*mu0hat/aI) &
- W_Y/larg/(1.+W_Y)*dlargdl * dmu0hat_aIdl &
+ (1.-W_Y)*d2mu0hat_aIdl2)
SI = 1. - max(exp(-aI*par*ThetaHat/mu0hat),0.)
dSIdl = ( (alphaI- dmu0hatdl/mu0hat) &
* ThetaHat + dThetaHatdl) * (1.-SI)*aI*par/mu0hat !confirmed
d2SIdl2 = par*(- dSIdl*aI/mu0hat*(ThetaHat*alphaI &
- ThetaHat/mu0hat*dmu0hatdl + dThetaHatdl) + (1.-SI) &
* (ThetaHat*alphaI- ThetaHat/mu0hat*dmu0hatdl + dThetaHatdl) &
* daI_mu0hatdl + (1.-SI)*aI/mu0hat*( &
- (d2mu0hatdl2/mu0hat - dmu0hatdl**2/mu0hat**2)*ThetaHat &
+ (alphaI-dmu0hatdl/mu0hat)*dThetaHatdl + d2ThetaHatdl2) )
! Light dependent growth rate
! (needs to take into account the cost of dark and light-dependent chl maintenance)
mu0hatSI = mu0hat*SI ! Gross specific carbon uptake (photosynthesis)
muIhat = mu0hatSI-(mu0hatSI+RMchl)*zetaChl*ThetaHat ! Net specific carbon uptake
muIhat = max(muIhat,eps)
dmu0hatSIdl = SI*dmu0hatdl + mu0hat*dSIdl
d2mu0hatSIdl2 = d2mu0hatdl2*SI+2.*dmu0hatdl*dSIdl+mu0hat*d2SIdl2 !Correct
dmuIhatdl = (1D0-zetaChl*ThetaHat)*dmu0hatSIdl - dThetaHatdl*zetaChl*(mu0hatSI+RMchl) !Correct
d2muIhatdl2=d2mu0hatSIdl2 - zetaChl*(ThetaHat*d2mu0hatSIdl2+2.*dThetaHatdl*dmu0hatSIdl &
+mu0hatSI*d2ThetaHatdl2)-zetaChl*RMchl*d2ThetaHatdl2 !Correct
ZINT = Qs*(muIhat/VNhat + zetaN)
dZINdl = Qs*(dmuIhatdl/VNhat - muIhat*dVNhatdl/VNhat**2)+alphaQ*ZINT
d2ZINdl2 = Qs/VNhat*((alphaQ-dVNhatdl/VNhat)*dmuIhatdl &
+ d2muIhatdl2) - Qs/VNhat**2*(muIhat*d2VNhatdl2 &
+ dVNhatdl*(dmuIhatdl+alphaQ*muIhat-2.*muIhat &
/ VNhat*dVNhatdl)) + alphaQ*dZINdl
fV = (-1.0 + sqrt(1.0 + 1.0/ZINT))*Qs*muIhat/VNhat
fV = max(fV,fVmin)
!
else
! Under the conditions of no light:
ThetaHat = thetamin ! a small positive value
dThetaHatdl = 0.
d2ThetaHatdl2 = 0.
ZINT = Qs*zetaN
dZINdl = alphaQ*ZINT
d2ZINdl2 = alphaQ*dZINdl
fV = fVmin
muIhat = 0.
dmuIhatdl = 0.
d2muIhatdl2 = 0.
endif
! Optimal nutrient quota:
QN = (1.+ sqrt(1.+1./ZINT))*Qs
dQNdl = alphaQ*QN-dZINdl*Qs/(2.*ZINT*sqrt(ZINT*(1.+ZINT))) !confirmed
d2QNdl2 = alphaQ*dQNdl - Qs/(2.*ZINT*sqrt(ZINT*(ZINT+1.))) &
*(d2ZINdl2+alphaQ*dZINdl-(2.*ZINT+1.5)/(ZINT*(ZINT+1.))*dZINdl**2) ! Confirmed
dfVdl = alphaQ*Qs*(1/QN+2.*zetaN)-(zetaN+Qs/QN**2)*dQNdl !Confirmed
!
d2fVdl2 = (alphaQ**2)*Qs*(1/QN + 2*zetaN) &
- 2.*alphaQ*Qs*dQNdl/QN**2 + 2.*(dQNdl**2)*Qs/QN**3 &
- (zetaN + Qs/QN**2) * d2QNdl2 ! Confirmed
if (par .gt. I_zero) then
!Net growth rate (d-1) of phytoplankton at the average size
muNet = muIhat*(1.-fV-Qs/QN) - zetaN*fV*VNhat
!Here the derivative of muNet includes respiratory costs of both N Assim and Chl maintenance
dmuNetdl = muIhat*(Qs/(QN**2)*dQNdl &
- alphaQ*Qs/QN-dfVdl) + (1.-fV-Qs/QN)*dmuIhatdl &
- zetaN*(fV*dVNhatdl+VNhat*dfVdl)
d2muNetdl2 = (Qs/(QN*QN))*dQNdl*dmuIhatdl &
+ muIhat*(Qs/QN**2)*(dQNdl*(alphaQ - 2.*dQNdl/QN) + d2QNdl2) &
- (alphaQ*(Qs/QN*(dmuIhatdl + alphaQ*muIhat) &
- muIhat*Qs/(QN**2)*dQNdl)) &
- (muIhat*d2fVdl2+dfVdl*dmuIhatdl) &
+ dmuIhatdl*(Qs/(QN**2)*dQNdl-alphaQ*Qs/QN-dfVdl) &
+ (1.-fV-Qs/QN)*d2muIhatdl2 &
- zetaN*(fV*d2VNhatdl2 + 2.*dfVdl*dVNhatdl + VNhat*d2fVdl2) !dC
else
muNet = 0.
dmuNetdl = 0.
d2muNetdl2 = 0.
endif
! chl:C ratio [ g chl / mol C ] of the whole cell, at the mean cell size
Theta = ThetaHat*(1. - fV - Qs/QN)
! Calculate the mean chl:C ratio of phytoplankton (averaged over all sizes)
! using the second derivative of the cell quota w.r.t. log size. (Why negative?)
ThetaAvg = max(Theta + VAR*.5*(d2ThetaHatdl2*(1.0 - fV - Qs/QN) &
- ThetaHat*(d2fVdl2 + Qs*(dQNdl**2*(2./QN)-d2QNdl2)/(QN**2))),thetamin)
! Calculate the mean N:C ratio of phytoplankton (averaged over all sizes)
! using the second derivative of the cell quota w.r.t. log size. (How to derive?)
QNAvg = max(QN/(1.+((2/QN)*dQNdl**2 - d2QNdl2)*VAR/(2.*QN)),QNmin)
!=============================================================
! Calculate the community sinking rate of phytoplankton
! Phytoplankton sinking rate at the average size
w_p = ScaleTrait(PMU,abs(dtdays*w_p0),alphaW) !Positive
! Constrain the sinking rate not too large for large cells (Smayda 1970)
w_p = min(w_p, Wmax*dtdays*1.)
dwdl = alphaW*w_p
d2wdl2 = alphaW*dwdl
! Sinking rate (m) of phytoplankton
w_pAvg = w_p+0.5*VAR*d2wdl2
! sinking rate must be negative!
w_pAvg = -min(w_pAvg, Wmax*dtdays*1.)
!=============================================================
!! ZOOplankton section:
tf_z = temp(Ez,tC)
! The grazing dependence on total prey
gbar = grazing(grazing_formulation,kp,PHY)
!Zooplankton total ingestion rate
INGES = tf_z*dtdays*gmax*gbar
!Zooplankton excretion rate (-> DOM)
RES = INGES*(1.-GGE-unass)
!ZOOPLANKTON EGESTION (-> POM)
EGES = INGES*unass
! Grazing rate on the mean size of PHY (specific to N-based Phy biomass, unit: d-1) (Eq. 12)
gbar = INGES*ZOO/PHY*sqrt(alphaG)
dgdlbar = 0.
d2gdl2bar = (1.-alphaG)/VAR*gbar
!!End of zooplankton section
!=============================================================
!! Solve ODE functions:
Zmort = ZOO*ZOO*dtdays*mz*tf_z !Mortality term for ZOO
if (d2gdl2bar .lt. 0.) then
VAR1 = VAR*(1.-VAR*d2gdl2bar)
else
VAR1 = VAR/(1.+VAR*d2gdl2bar)
endif
! Update PMU and VAR:
if (d2muNetdl2 .gt. 0.) then
! self%dVARdt=VAR*VAR*(self%d2muNetdl2- self%d2gdl2bar-self%d2wdl2) !Eq. 19
VAR1 =VAR1*(1.+VAR*d2muNetdl2)
else
VAR1 =VAR1/(1.-VAR*d2muNetdl2)
endif
VAR1 =VAR1/(1.+VAR*d2wdl2)
!Eq. 18, Update PMU:
! self%dPMUdt = self%p%VAR*(self%dmuNetdl - self%dgdlbar - self%dwdl)
if (dmuNetdl .gt. 0.) then
PMU1 = PMU + VAR * dmuNetdl
else
PMU1 = PMU/(1.-VAR/PMU*dmuNetdl)
endif
PMU1 = PMU1/(1.+ VAR/PMU*dwdl)
!Contrain the PMU and VAR:
PMU1 = min(max(PMU1,0.01),PMUmax)
VAR1 = min(max(VAR1,0.01),VARmax)
! For production/destruction matrix:
pp(:,:) = 0.
pp(iNO3,iDET)=dtdays*rDN*DET*tf_z
pp(iNO3,iZOO)=ZOO*RES
pp(iDET,iZOO)=ZOO*EGES+Zmort
pp(iZOO,iPHY)=ZOO*INGES
pp(iDET,iPHY)=PHY*PHY*tf_p*dtdays*mp
PP_pn = PHY*(muNet + 0.5*VAR*d2muNetdl2)
pp(iPHY,iNO3)=max(0.5*( PP_pn + abs(PP_pn)), 0.)
pp(iNO3,iPHY)=max(0.5*(-PP_pn + abs(PP_pn)), 0.)
DET1 = (DET + pp(iDET,iZOO) + pp(iDET,iPHY))/(1.+ pp(iNO3,iDET)/DET )
DET1 = min(max(DET1,eps),DETmax)
NO31 = (NO3+pp(iNO3,iDET)+pp(iNO3,iZOO)+pp(iNO3,iPHY))/(1.+pp(iPHY,iNO3)/NO3)
PHY1=(PHY+pp(iPHY,iNO3)) &
/(1.+(pp(iZOO,iPHY)+pp(iNO3,iPHY)+pp(iDET,iPHY))/PHY)
ZOO1 = (ZOO+pp(iZOO,iPHY))/(1.+ EGES+ZOO*dtdays*mz*tf_z+RES)
CHL1 = PHY1/QNavg*ThetaAvg
PMUPHY1 =(PMUPHY+PHY*PMU1+PMU*PHY1)/(1.+ 2.*PHY*PMU/PMUPHY)
VARPHY1 =(VARPHY+PHY*VAR1+VAR*PHY1)/(1.+ 2.*PHY*VAR/VARPHY)
Varout(oNO3) = NO31
Varout(oPHY) = PHY1
Varout(oZOO) = ZOO1
Varout(oDET) = DET1
Varout(oPMU) = PMUPHY1
Varout(oVAR) = VARPHY1
Varout(oCHL) = CHL1
if(do_IRON .eq. .true.) Varout(oFER) = Fe
Varout(oTheta) = ThetaAvg
Varout(oQN) = QNAvg
Varout(omuNet) = PP(iPHY,iNO3)/dtdays/PHY
Varout(oGraz ) = PP(iZOO,iPHY)/dtdays/PHY
Varout(oZ2N ) = PP(iNO3,iZOO)/dtdays
Varout(oD2N ) = PP(iNO3,iDET)/dtdays
Varout(odmudl) = dmuNetdl/dtdays
Varout(od2mu ) = d2muNetdl2/dtdays
Varout(od2gdl) = d2gdl2bar/dtdays
Varout(ow_p) = w_pAvg/dtdays
return
END SUBROUTINE FlexEFT
!========================================================
real(8) function temp(Ea,tC)
implicit none
!DESCRIPTION:
!The temperature dependence of plankton rates are fomulated according to the Arrhenuis equation.
! tC: in situ temperature
! Tr: reference temperature
!
!INPUT PARAMETERS:
real(8), intent (in) :: Ea, tC
! boltzman constant constant [ eV /K ]
real(8), parameter :: kb = 8.62d-5, Tr = 15.
temp = exp(-(Ea/kb)*(1D0/(273.15 + tC)-1D0/(273.15 + Tr)))
return
end function temp
!========================================================
FUNCTION WAPR (X, NB, L) RESULT (WAP)
!
! WAPR - output
! X - argument of W(X)
! NB is the branch of the W function needed:
! NB = 0 - upper branch
! NB <> 0 - lower branch
!
! NERROR is the output error flag:
! NERROR = 0 -> routine completed successfully
! NERROR = 1 -> X is out of range
!
! Range: -exp(-1) <= X for the upper branch of the W function
! -exp(-1) < X < 0 for the lower branch of the W function
!
! L - determines how WAPR is to treat the argument X
! L = 1 -> X is the offset from -exp(-1), so compute
! W(X-exp(-1))
! L <> 1 -> X is the desired X, so compute W(X)
!
! M - print messages from WAPR?
! M = 1 -> Yes
! M <> 1 -> No
!
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! NN is the output device number
!
! NBITS is the number of bits (less 1) in the mantissa of the
! floating point number number representation of your machine.
! It is used to determine the level of accuracy to which the W
! function should be calculated.
!
! Most machines use a 24-bit matissa for single precision and
! 53-56 bits for double precision. The IEEE standard is 53
! bits. The Fujitsu VP2200 uses 56 bits. Long word length
! machines vary, e.g., the Cray X/MP has a 48-bit mantissa for
! single precision.
!
IMPLICIT NONE
real, INTENT(in) :: X
INTEGER, PARAMETER :: NN=6, NBITS=23, NITER=1
real, PARAMETER ::EM=-0.367879441171442,& ! -EXP(-1)
EM9=-1.234098040866796E-4,& ! -EXP(-9)
C13=1.E0/3.E0,&
C23=2.E0*C13,&
EM2=2.E0/EM,&
E12=-EM2,&
TB=.5E0**NBITS,&
TB2=.5E0**(NBITS/2),& ! SQRT(TB)
X0=0.0350769390096679055,& ! TB**(1/6)*0.5E0
X1=0.302120119432784731,& !(1 - 17*TB**(2/7))*EM
AN22=3.6E2/83.E0,&
AN11=135./83.E0,&
AN3=8.E0/3.E0,&
AN4=135.E0/83.E0,&
AN5=166.E0/39.E0,&
AN6=3167.E0/3549.E0,&
S2=1.41421356237310,& ! SQRT(2.E0)
S21=2.E0*S2-3.E0,&
S22=4.E0-3.E0*S2,&
S23=S2-2.E0
real :: WAP, AN2, DELX, RETA, ZL, T, TS, ETA, TEMP, TEMP2, ZN
real :: XX
INTEGER, INTENT(in) :: NB, L
! Various mathematical constants
!
! The following COMMON statement is needed only when testing this
! function using BISECT, otherwise it can be removed.
!
! COMMON/WAPCOM/NBITS
! DATA INIT,NITER/0,1/
! DATA NBITS/23/
!
! IF(INIT.EQ.0) THEN
! INIT=1
!
! Code to calculate NBITS for the host machine. NBITS is the
! mantissa length less one. This value is chosen to allow for
! rounding error in the final bit. This need only be run once on
! any particular machine. It can then be included in the above
! DATA statement.
!
! DO I=1,2000
! B=2.E0**(-I)
! V=1.E0+B
! IF(V.EQ.1.E0)THEN
! NBITS=I-1
! J=-ALOG10(B)
! IF(M.EQ.1) WRITE(NN,40)NBITS,J
! EXIT
! ENDIF
! END DO