Commit 8ea84af1b012817d8e1463ca5b09810f6ba5e354

Authored by dumoda01
1 parent e74c0f0a
Exists in master and in 1 other branch snow

Correction d'un bogue dans bio.F90 pour que le code fonctionne avec bio_ismer.F9…

…0. Pour changer de modele biologique, il faut modifier l'indice des traceurs nit et amm dans bio.F90 et meanflow/nitrate.F90 et meanflow/ammonium.F90. Ceci devra etre regularise eventuellement.
Showing 2 changed files with 128 additions and 133 deletions   Show diff stats
src/extras/bio/bio.F90
@@ -565,11 +565,11 @@ @@ -565,11 +565,11 @@
565 if (bio_eulerian) then 565 if (bio_eulerian) then
566 do j=1,numcc 566 do j=1,numcc
567 567
568 - if(j==5) then !CHG3 568 + if(j==1) then !CHG3
569 do k=1,nlev 569 do k=1,nlev
570 cc(j,k) = nit(k) 570 cc(j,k) = nit(k)
571 end do 571 end do
572 - else if (j==6) then !CHG5 572 + else if (j==9) then !CHG5
573 do k=1,nlev 573 do k=1,nlev
574 cc(j,k) = amm(k) 574 cc(j,k) = amm(k)
575 end do 575 end do
src/extras/bio/bio_ismer.F90
1 -!$Id: bio_fasham.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $ 1 +!$Id: bio_ismer.F90,v 1.11 2007-01-06 11:49:15 kbk Exp $
2 #include"cppdefs.h" 2 #include"cppdefs.h"
3 !----------------------------------------------------------------------- 3 !-----------------------------------------------------------------------
4 !BOP 4 !BOP
@@ -54,61 +54,58 @@ @@ -54,61 +54,58 @@
54 ! 54 !
55 ! !LOCAL VARIABLES: 55 ! !LOCAL VARIABLES:
56 ! from a namelist 56 ! from a namelist
57 - REALTYPE :: p1_initial= 0.05  
58 - REALTYPE :: p2_initial= 0.05  
59 - REALTYPE :: z1_initial= 0.05  
60 - REALTYPE :: z2_initial= 0.05  
61 - REALTYPE :: b_initial= 0.001  
62 - REALTYPE :: d_initial= 0.416666666  
63 -! Nitrate and ammonium are initialized within the GOTM observation module  
64 -! REALTYPE :: n_initial= 8.3  
65 -! REALTYPE :: a_initial= 0.22  
66 - REALTYPE :: l_initial= 0.14  
67 - REALTYPE :: p0 = 0.0  
68 - REALTYPE :: z0 = 0.0  
69 - REALTYPE :: b0 = 0.0  
70 - REALTYPE :: vp1 = 1.5  
71 - REALTYPE :: alpha1 = 0.065  
72 - REALTYPE :: inib1 = 0.05  
73 - REALTYPE :: vp2 = 1.5  
74 - REALTYPE :: alpha2 = 0.065  
75 - REALTYPE :: inib2 = 0.05  
76 - REALTYPE :: theta = 0.0  
77 - REALTYPE :: w_p1min = -0.06  
78 - REALTYPE :: w_p1max = -0.38  
79 - REALTYPE :: w_p2min = -0.06  
80 - REALTYPE :: w_p2max = -0.38  
81 - REALTYPE :: kn1 = 0.2  
82 - REALTYPE :: ka1 = 0.8  
83 - REALTYPE :: kn2 = 0.2  
84 - REALTYPE :: ka2 = 0.8  
85 - REALTYPE :: mu1 = 0.05  
86 - REALTYPE :: k5 = 0.2  
87 - REALTYPE :: gamma = 0.05  
88 - REALTYPE :: w_p1 = -0.5  
89 - REALTYPE :: w_p2 = -0.5  
90 - REALTYPE :: g1max = 1.0  
91 - REALTYPE :: g2max = 1.0  
92 - REALTYPE :: k3 = 1.0  
93 - REALTYPE :: beta = 0.625  
94 - REALTYPE :: mu2 = 0.3  
95 - REALTYPE :: k6 = 0.2  
96 - REALTYPE :: delta = 0.1  
97 - REALTYPE :: epsi = 0.70  
98 - REALTYPE :: r11 = 0.55  
99 - REALTYPE :: r12 = 0.4  
100 - REALTYPE :: r13 = 0.05  
101 - REALTYPE :: r21 = 0.55  
102 - REALTYPE :: r22 = 0.4  
103 - REALTYPE :: r23 = 0.05  
104 - REALTYPE :: vb = 1.2  
105 - REALTYPE :: k4 = 0.5  
106 - REALTYPE :: mu3 = 0.15  
107 - REALTYPE :: eta = 0.0  
108 - REALTYPE :: mu4 = 0.02  
109 - REALTYPE :: mu5 = 0.00  
110 - REALTYPE :: w_d = -2.0  
111 - REALTYPE, public :: kc = 0.03 57 + REALTYPE :: p1_init = 0.05
  58 + REALTYPE :: p2_init = 0.05
  59 + REALTYPE :: z1_init = 0.05
  60 + REALTYPE :: z2_init = 0.05
  61 + REALTYPE :: b_init = 0.001
  62 + REALTYPE :: d_init = 0.4
  63 + REALTYPE :: l_init = 0.14
  64 + REALTYPE :: p0 = 0.0
  65 + REALTYPE :: z0 = 0.0
  66 + REALTYPE :: b0 = 0.0
  67 + REALTYPE :: vp1 = 1.5
  68 + REALTYPE :: alpha1 = 0.065
  69 + REALTYPE :: inib1 = 0.05
  70 + REALTYPE :: vp2 = 1.5
  71 + REALTYPE :: alpha2 = 0.065
  72 + REALTYPE :: inib2 = 0.05
  73 + REALTYPE :: theta = 0.0
  74 + REALTYPE :: w_p1min = -0.06
  75 + REALTYPE :: w_p1max = -0.38
  76 + REALTYPE :: w_p2min = -0.06
  77 + REALTYPE :: w_p2max = -0.38
  78 + REALTYPE :: kn1 = 0.2
  79 + REALTYPE :: ka1 = 0.8
  80 + REALTYPE :: kn2 = 0.2
  81 + REALTYPE :: ka2 = 0.8
  82 + REALTYPE :: mu1 = 0.05
  83 + REALTYPE :: k5 = 0.2
  84 + REALTYPE :: gamma = 0.05
  85 + REALTYPE :: w_p1 = -0.5
  86 + REALTYPE :: w_p2 = -0.5
  87 + REALTYPE :: g1max = 1.0
  88 + REALTYPE :: g2max = 1.0
  89 + REALTYPE :: k3 = 1.0
  90 + REALTYPE :: beta = 0.625
  91 + REALTYPE :: mu2 = 0.3
  92 + REALTYPE :: k6 = 0.2
  93 + REALTYPE :: delta = 0.1
  94 + REALTYPE :: epsi = 0.70
  95 + REALTYPE :: r11 = 0.55
  96 + REALTYPE :: r12 = 0.4
  97 + REALTYPE :: r13 = 0.05
  98 + REALTYPE :: r21 = 0.55
  99 + REALTYPE :: r22 = 0.4
  100 + REALTYPE :: r23 = 0.05
  101 + REALTYPE :: vb = 1.2
  102 + REALTYPE :: k4 = 0.5
  103 + REALTYPE :: mu3 = 0.15
  104 + REALTYPE :: eta = 0.0
  105 + REALTYPE :: mu4 = 0.02
  106 + REALTYPE :: mu5 = 0.00
  107 + REALTYPE :: w_d = -2.0
  108 + REALTYPE, public :: kc = 0.03
112 integer :: out_unit 109 integer :: out_unit
113 integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b=8,a=9 110 integer, parameter :: n=1,p1=2,p2=3,z1=4,z2=5,d=6,l=7,b=8,a=9
114 !EOP 111 !EOP
@@ -143,14 +140,14 @@ @@ -143,14 +140,14 @@
143 ! 140 !
144 ! !LOCAL VARIABLES: 141 ! !LOCAL VARIABLES:
145 namelist /bio_ismer_nml/ numc, & 142 namelist /bio_ismer_nml/ numc, &
146 - p1_initial,p2_initial,z1_initial,z2_initial, &  
147 - b_initial,d_initial,l_initial, &  
148 - p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &  
149 - kn1,ka1,kn2,ka2,mu1,k5,gamma,w_p1,w_p2, &  
150 - g1max,g2max,k3,beta,mu2,k6,delta,epsi, &  
151 - r11,r12,r13,r21,r22,r23, &  
152 - vb,k4,mu3,eta,mu4,w_d,kc,mu5, &  
153 - theta,w_p1max,w_p1min,w_p2min,w_p2max !CHG2 143 + p1_init,p2_init,z1_init,z2_init, &
  144 + b_init,d_init,l_init, &
  145 + p0,z0,b0,vp1,alpha1,inib1,vp2,alpha2,inib2, &
  146 + kn1,ka1,kn2,ka2,mu1,k5,gamma,w_p1,w_p2, &
  147 + g1max,g2max,k3,beta,mu2,k6,delta,epsi, &
  148 + r11,r12,r13,r21,r22,r23, &
  149 + vb,k4,mu3,eta,mu4,w_d,kc,mu5, &
  150 + theta,w_p1max,w_p1min,w_p2min,w_p2max
154 !EOP 151 !EOP
155 !----------------------------------------------------------------------- 152 !-----------------------------------------------------------------------
156 !BOC 153 !BOC
@@ -178,28 +175,28 @@ @@ -178,28 +175,28 @@
178 901 format (f8.5) 175 901 format (f8.5)
179 176
180 ! Conversion from day to second 177 ! Conversion from day to second
181 - vp1 = vp1 /secs_pr_day  
182 - vp2 = vp2 /secs_pr_day  
183 - vb = vb /secs_pr_day  
184 - mu1 = mu1 /secs_pr_day  
185 - mu2 = mu2 /secs_pr_day  
186 - mu3 = mu3 /secs_pr_day  
187 - mu4 = mu4 /secs_pr_day  
188 - mu5 = mu5 /secs_pr_day !DD  
189 - g1max = g1max /secs_pr_day  
190 - g2max = g2max /secs_pr_day  
191 - w_p1 = w_p1 /secs_pr_day  
192 - w_p2 = w_p2 /secs_pr_day  
193 - w_p1min = w_p1min /secs_pr_day !DD  
194 - w_p1max = w_p1max /secs_pr_day !DD  
195 - w_p2min = w_p2min /secs_pr_day !DD  
196 - w_p2max = w_p2max /secs_pr_day !DD  
197 - theta = theta /secs_pr_day !DD  
198 - w_d = w_d /secs_pr_day  
199 - alpha1 = alpha1/secs_pr_day  
200 - inib1 = inib1 /secs_pr_day !CHG1  
201 - alpha2 = alpha2/secs_pr_day  
202 - inib2 = inib2 /secs_pr_day !CHG1 178 + vp1 = vp1 /secs_pr_day
  179 + vp2 = vp2 /secs_pr_day
  180 + vb = vb /secs_pr_day
  181 + mu1 = mu1 /secs_pr_day
  182 + mu2 = mu2 /secs_pr_day
  183 + mu3 = mu3 /secs_pr_day
  184 + mu4 = mu4 /secs_pr_day
  185 + mu5 = mu5 /secs_pr_day
  186 + g1max = g1max /secs_pr_day
  187 + g2max = g2max /secs_pr_day
  188 + w_p1 = w_p1 /secs_pr_day
  189 + w_p2 = w_p2 /secs_pr_day
  190 + w_p1min = w_p1min /secs_pr_day
  191 + w_p1max = w_p1max /secs_pr_day
  192 + w_p2min = w_p2min /secs_pr_day
  193 + w_p2max = w_p2max /secs_pr_day
  194 + theta = theta /secs_pr_day
  195 + w_d = w_d /secs_pr_day
  196 + alpha1 = alpha1 /secs_pr_day
  197 + inib1 = inib1 /secs_pr_day
  198 + alpha2 = alpha2 /secs_pr_day
  199 + inib2 = inib2 /secs_pr_day
203 200
204 out_unit=unit 201 out_unit=unit
205 202
@@ -250,13 +247,13 @@ @@ -250,13 +247,13 @@
250 !BOC 247 !BOC
251 do i=1,nlev 248 do i=1,nlev
252 cc(n,i) = nprof(i) !CHG3 249 cc(n,i) = nprof(i) !CHG3
253 - cc(p1,i)= p1_initial  
254 - cc(p2,i)= p2_initial  
255 - cc(z1,i)= z1_initial  
256 - cc(z2,i)= z2_initial  
257 - cc(d,i) = d_initial  
258 - cc(l,i) = l_initial  
259 - cc(b,i) = b_initial 250 + cc(p1,i)= p1_init
  251 + cc(p2,i)= p2_init
  252 + cc(z1,i)= z1_init
  253 + cc(z2,i)= z2_init
  254 + cc(d,i) = d_init
  255 + cc(l,i) = l_init
  256 + cc(b,i) = b_init
260 cc(a,i) = aprof(i) !CHG5 257 cc(a,i) = aprof(i) !CHG5
261 end do 258 end do
262 259
@@ -506,7 +503,7 @@ @@ -506,7 +503,7 @@
506 ! 503 !
507 ! !LOCAL VARIABLES: 504 ! !LOCAL VARIABLES:
508 REALTYPE :: fac1,fac2,fac3,minal,qn1,qa1,qn2,qa2 505 REALTYPE :: fac1,fac2,fac3,minal,qn1,qa1,qn2,qa2
509 - REALTYPE :: Ps1,Ps2,ff1,ff2 !CHG1 506 + REALTYPE :: ps1,ps2,ff1,ff2 !CHG1
510 integer :: i,j,ci 507 integer :: i,j,ci
511 !EOP 508 !EOP
512 !----------------------------------------------------------------------- 509 !-----------------------------------------------------------------------
@@ -531,19 +528,22 @@ @@ -531,19 +528,22 @@
531 ! Parker (1974) - inhibition 528 ! Parker (1974) - inhibition
532 ! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))**2 529 ! ff= vp*((par(ci)/I_opt)*exp(1-(par(ci)/I_opt)))**2
533 ! Platt et al. (1980) - inhibition 530 ! Platt et al. (1980) - inhibition
534 - Ps1= vp1/((alpha1/(alpha1+inib1))*(alpha1/(alpha1+inib1))**(inib1/alpha1))  
535 - ff1= Ps1*(1.-exp(-1.*alpha1*par(ci)/Ps1))*exp(-1.*inib1*par(ci)/Ps1)  
536 - Ps2= vp2/((alpha2/(alpha2+inib2))*(alpha2/(alpha2+inib2))**(inib2/alpha2))  
537 - ff2= Ps2*(1.-exp(-1.*alpha2*par(ci)/Ps2))*exp(-1.*inib2*par(ci)/Ps2)  
538 -! --------------------------------------------------------------------  
539 531
  532 + ! Light limitation factor (PI curve)
  533 + ps1= vp1/((alpha1/(alpha1+inib1))*(alpha1/(alpha1+inib1))**(inib1/alpha1))
  534 + ff1= ps1*(1.-exp(-1.*alpha1*par(ci)/ps1))*exp(-1.*inib1*par(ci)/ps1)
540 535
  536 + ps2= vp2/((alpha2/(alpha2+inib2))*(alpha2/(alpha2+inib2))**(inib2/alpha2))
  537 + ff2= ps2*(1.-exp(-1.*alpha2*par(ci)/ps2))*exp(-1.*inib2*par(ci)/ps2)
  538 +
  539 + ! Nutrient limitation factors
541 qn1=(cc(n,ci)/kn1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1) 540 qn1=(cc(n,ci)/kn1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
542 qa1=(cc(a,ci)/ka1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1) 541 qa1=(cc(a,ci)/ka1)/(1.+cc(n,ci)/kn1+cc(a,ci)/ka1)
543 542
544 qn2=(cc(n,ci)/kn2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2) 543 qn2=(cc(n,ci)/kn2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
545 qa2=(cc(a,ci)/ka2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2) 544 qa2=(cc(a,ci)/ka2)/(1.+cc(n,ci)/kn2+cc(a,ci)/ka2)
546 545
  546 + ! Grazing preference normalization factors
547 fac1=(cc(z1,ci)+z0)/(k3*(r11*cc(p1,ci)+r12*cc(b,ci)+r13*cc(d,ci))+ & 547 fac1=(cc(z1,ci)+z0)/(k3*(r11*cc(p1,ci)+r12*cc(b,ci)+r13*cc(d,ci))+ &
548 r11*cc(p1,ci)**2+r12*cc(b,ci)**2+r13*cc(d,ci)**2) 548 r11*cc(p1,ci)**2+r12*cc(b,ci)**2+r13*cc(d,ci)**2)
549 549
@@ -562,49 +562,44 @@ @@ -562,49 +562,44 @@
562 nitlim(ci) =qn1 562 nitlim(ci) =qn1
563 ammlim(ci) =qa1 563 ammlim(ci) =qa1
564 564
565 - dd(p1,d,ci)=mu1*(cc(p1,ci)+p0)/(k5+cc(p1,ci)+p0)*cc(p1,ci) &  
566 - +fac3*(1.-beta)*cc(p1,ci)**2*(g1max*r11*fac1 + g2max*r21*fac2)  
567 - dd(p2,d,ci)=mu1*(cc(p2,ci)+p0)/(k5+cc(p2,ci)+p0)*cc(p2,ci) &  
568 - +fac3*(1.-beta)*g2max*r12*cc(p2,ci)**2*fac2 565 + ! Nutrient uptake by flagellates and diatoms
  566 + dd(n,p1,ci) =ff1*qn1*(cc(p1,ci)+p0)
  567 + dd(a,p1,ci) =ff1*qa1*(cc(p1,ci)+p0)
  568 + dd(n,p2,ci) =ff2*qn2*(cc(p2,ci)+p0)
  569 + dd(a,p2,ci) =ff2*qa2*(cc(p2,ci)+p0)
  570 +
  571 + dd(p1,l,ci) =gamma*ff1*(qn1+qa1)*cc(p1,ci)
  572 + dd(p2,l,ci) =gamma*ff2*(qn2+qa2)*cc(p2,ci)
569 573
570 - dd(p1,l,ci)=gamma*ff1*(qn1+qa1)*cc(p1,ci)  
571 - dd(p2,l,ci)=gamma*ff2*(qn2+qa2)*cc(p2,ci) 574 + dd(p1,d,ci) =mu1*(cc(p1,ci)+p0)/(k5+cc(p1,ci)+p0)*cc(p1,ci) &
  575 + +fac3*(1.-beta)*cc(p1,ci)**2*(g1max*r11*fac1+g2max*r21*fac2)
  576 + dd(p2,d,ci) =mu1*(cc(p2,ci)+p0)/(k5+cc(p2,ci)+p0)*cc(p2,ci) &
  577 + +fac3*(1.-beta)*g2max*r21*cc(p2,ci)**2*fac2
572 578
573 - dd(b,d,ci) =fac3*(1.-beta)*g1max*r12*cc(b,ci)**2*fac1 579 + dd(b,d,ci) =fac3*(1.-beta)*g1max*r12*cc(b,ci)**2*fac1
574 580
575 dd(p1,z1,ci)=fac3*beta*g1max*r11*cc(p1,ci)**2*fac1 581 dd(p1,z1,ci)=fac3*beta*g1max*r11*cc(p1,ci)**2*fac1
576 - dd(b,z1,ci)=fac3*beta*g1max*r12*cc(b,ci)**2*fac1  
577 - dd(d,z1,ci)=fac3*beta*g1max*r13*cc(d,ci)**2*fac1 582 + dd(b,z1,ci) =fac3*beta*g1max*r12*cc(b,ci)**2*fac1
  583 + dd(d,z1,ci) =fac3*beta*g1max*r13*cc(d,ci)**2*fac1
578 584
579 dd(p1,z2,ci)=fac3*beta*g2max*r21*cc(p1,ci)**2*fac2 585 dd(p1,z2,ci)=fac3*beta*g2max*r21*cc(p1,ci)**2*fac2
580 dd(p2,z2,ci)=fac3*beta*g2max*r22*cc(p2,ci)**2*fac2 586 dd(p2,z2,ci)=fac3*beta*g2max*r22*cc(p2,ci)**2*fac2
581 - dd(d,z2,ci)=fac3*beta*g2max*r23*cc(d,ci)**2*fac2 587 + dd(d,z2,ci) =fac3*beta*g2max*r23*cc(d,ci)**2*fac2
582 588
583 - dd(b,a,ci) =mu3*cc(b,ci)  
584 - dd(d,l,ci) =mu4*cc(d,ci)  
585 - dd(a,n,ci) =mu5*cc(a,ci) 589 + dd(b,a,ci) =mu3*cc(b,ci)
  590 + dd(d,l,ci) =mu4*cc(d,ci)
  591 + dd(a,n,ci) =mu5*cc(a,ci)
586 592
587 - dd(z1,d,ci)=(1.-epsi-delta)*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)  
588 - dd(z1,a,ci)=epsi*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)  
589 - dd(z1,l,ci)=delta*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)  
590 -  
591 - dd(z2,d,ci)=(1.-epsi-delta)*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)  
592 - dd(z2,a,ci)=epsi*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)  
593 - dd(z2,l,ci)=delta*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)  
594 -  
595 - dd(n,p1,ci) =ff1*qn1*(cc(p1,ci)+p0)  
596 - dd(a,p1,ci) =ff1*qa1*(cc(p1,ci)+p0)  
597 - dd(n,p2,ci) =ff2*qn2*(cc(p2,ci)+p0)  
598 - dd(a,p2,ci) =ff2*qa2*(cc(p2,ci)+p0) 593 + dd(z1,d,ci) =(1.-epsi-delta)*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  594 + dd(z1,a,ci) =epsi*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
  595 + dd(z1,l,ci) =delta*mu2*(cc(z1,ci)+z0)/(k6+cc(z1,ci)+z0)*cc(z1,ci)
599 596
600 - dd(a,b,ci) =vb*minal/(k4+minal+cc(l,ci))*(cc(b,ci)+b0)  
601 - dd(l,b,ci) =vb*cc(l,ci)/(k4+minal+cc(l,ci))*(cc(b,ci)+b0) 597 + dd(z2,d,ci) =(1.-epsi-delta)*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  598 + dd(z2,a,ci) =epsi*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
  599 + dd(z2,l,ci) =delta*mu2*(cc(z2,ci)+z0)/(k6+cc(z2,ci)+z0)*cc(z2,ci)
602 600
603 -! Taux de chute du phytoplancton en fonction de la croissance  
604 -! ws(p,ci)=w_pmin+(w_pmax-w_pmin)*exp(-ff/theta)  
605 -! ws(p,ci)=w_pmax-(w_pmax-w_pmin)*(ff/theta/(1.+ff/theta))  
606 - ws(p1,ci)=w_p1  
607 - ws(p2,ci)=w_p2 601 + dd(a,b,ci) =vb*minal/(k4+minal+cc(l,ci))*(cc(b,ci)+b0)
  602 + dd(l,b,ci) =vb*cc(l,ci)/(k4+minal+cc(l,ci))*(cc(b,ci)+b0)
608 603
609 do i=1,numc 604 do i=1,numc
610 do j=1,numc 605 do j=1,numc