Index: NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90	(revision 10222)
+++ NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90	(revision 10362)
@@ -591,6 +591,10 @@
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
             zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4)
-            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
             zsedtra(5) = zbeta  - zsedtra(4)
@@ -601,6 +605,10 @@
             zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
-            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(2)
             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2)
@@ -609,6 +617,10 @@
             zalpha = zsedtra(1) - 2.0 * zsedtra(3)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
-            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(3)
             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
@@ -617,6 +629,10 @@
             zbeta  = zsedtra(4) + zsedtra(6)
             zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
-            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(6)
             zsedtra(4) = zbeta  - zsedtra(6)
@@ -626,6 +642,10 @@
             zbeta  = zsedtra(4) + zsedtra(6)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
-            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(2) = zalpha + zsedtra(4)
             zsedtra(6) = zbeta  - zsedtra(4)
@@ -637,6 +657,10 @@
             zdelta = pwcp(ji,jk,jwso4) + zsedtra(2)
             zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5)
-            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  &
-            &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 )
+            ELSE
+               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  &
+               &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) )
+            ENDIF
             zsedtra(5) = zalpha + 2.0 * zsedtra(2)
             zsedtra(4) = zbeta  - zsedtra(5)
@@ -648,6 +672,10 @@
             zbeta  = zsedtra(4) + zsedtra(6)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
-            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(2) = zalpha + zsedtra(4)
             zsedtra(6) = zbeta  - zsedtra(4)
@@ -657,6 +685,10 @@
             zbeta  = zsedtra(4) + zsedtra(6)
             zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
-            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
+            IF (zalpha == 0.) THEN
+               zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. )
+            ELSE
+               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(6)
             zsedtra(4) = zbeta  - zsedtra(6)
@@ -665,6 +697,10 @@
             zalpha = zsedtra(1) - 2.0 * zsedtra(3)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
-            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
+            IF (zalpha == 0.) THEN
+               zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0) 
+            ELSE
+               zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(3)
             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
@@ -673,6 +709,10 @@
             zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
-            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 2.0 * zsedtra(2)
             pwcp(ji,jk,jwso4) = zbeta - zsedtra(2)
@@ -683,6 +723,10 @@
             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
             zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4)
-            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
-            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
+            IF ( zalpha == 0. ) THEN
+               zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 )
+            ELSE
+               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) )
+            ENDIF
             zsedtra(1) = zalpha + 0.25 * zsedtra(4)
             zsedtra(5) = zbeta  - zsedtra(4)
