source: trunk/SOURCES/source_3.20/ice_ikaite.f @ 4

Last change on this file since 4 was 4, checked in by vancop, 8 years ago

initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 6.3 KB
Line 
1      SUBROUTINE ice_ikaite(nlay_i)
2
3!------------------------------------------------------------------------------!
4!                               *** ice_ikaite ***
5!
6!     This routine computes precipitation / dissolution of ikaite
7!     and impacts on DIC, Alk and Ca2+
8!
9!     Ref: Moreau et al (JGR2015)
10!     
11!     Original code: S. Moreau, M. Vancoppenolle, 2012-2015
12!     Rewriting: M. Vancopppenolle, Nov 2015
13!
14!     Solubility product formulation is valid until -9C
15!     The fit is not monotonic below -9C, so we may experience problems
16!     with ikaite below -9C
17!
18!------------------------------------------------------------------------------!
19
20      INCLUDE 'type.com'
21      INCLUDE 'para.com'
22      INCLUDE 'const.com'
23      INCLUDE 'ice.com'
24      INCLUDE 'thermo.com'
25      INCLUDE 'bio.com'
26     
27       REAL(8), DIMENSION(nlay_bio) :: 
28     &    Kspi
29
30       LOGICAL ::
31     &   ln_write_bio = .TRUE.
32
33!==============================================================================!
34
35!
36!---------------------------------------------------------------------
37! 1) Starting the routine
38!---------------------------------------------------------------------
39!
40
41      IF ( ln_write_bio ) THEN
42
43         WRITE(numout,*)
44         WRITE(numout,*) ' *** ice_ikaite : '
45         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~'
46         WRITE(numout,*)
47         WRITE(numout,*) ' *** Initial values before CaCO3 precip/diss '
48         WRITE(numout,*) '        at time step : ', numit
49
50      ENDIF
51     
52      CALL ice_brine     ! get brine salinity and density
53     
54      CALL ice_carb_chem ! get CO32m
55
56      DO jn = 1, ntra_bio
57
58         ! Brine / Total gas concentrations
59         DO layer = 1, nlay_bio
60             c_i_bio(jn,layer)   = cbu_i_bio(jn,layer) / e_i_bio(layer)
61         END DO
62
63         IF ( ln_write_bio ) THEN
64
65            WRITE(numout,*)
66            WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
67            WRITE(numout,*) ' c_i_bio: ', ( c_i_bio(jn,layer),  ! concentration of tracers in bulk ice
68     &                                    layer = 1, nlay_bio )
69            WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer),  ! concentration of tracers in bulk ice
70     &                                      layer = 1, nlay_bio )
71            WRITE(numout,*) 
72         ENDIF
73      END DO
74     
75!
76!------------------------------------------------------------------------------!
77! 2) Solubility product
78!------------------------------------------------------------------------------!
79!
80      ! Ref: Papadimitriou et al. (GCA 2013)
81      ! fit valid until -9 degC
82      ! non monotonic (at all) below this value
83      ! we may have to regularize this
84     
85      DO layer = 1, nlay_bio
86
87          zT  = t_i_bio(layer)
88
89          IF ( zt .LE. 273.15-9.) zT = 273.15-9. ! asked stathys papadim what to do here
90
91          zpKsp = - 15489.09608 + 623443.70216 / zT +
92     &             2355.14596 * log(zT)
93
94          Kspi(layer) = 10.**(-zpKsp)     ! (mol/kg)^2 of brine
95         
96          ! Conversion into (mmol/m3)^2
97          Kspi(layer) = Kspi(layer) * 1.0e6 * 
98     &                  rhobr_bio(layer) * rhobr_bio(layer)
99
100      END DO
101       
102      WRITE(numout,*) ' Kspi : ', Kspi(1:nlay_bio)
103     
104!
105!------------------------------------------------------------------------------!
106! 3) Calcification / Dissolution
107!------------------------------------------------------------------------------!
108!
109
110      DO layer = 1, nlay_bio
111
112         zcal = c_i_bio(jn_cal,layer)
113         ! zcal = 311.765 * sbr_bio(layer) --> approach not valid because Ca2+ changes a lot in brine in winter
114
115         ! Brine concentration of CO32-
116         zco32m_br  = co32m(layer) / e_i_bio(layer)
117         zco32m_sat = Kspi(layer)  / zcal
118         
119         ! Saturation state (omega) in brine
120         zomega = zco32m_br * zcal / Kspi(layer)
121         
122         ! Precipitation / dissolution
123         zd_ika = ddtb / caco3_time * ( zco32m_br - zco32m_sat )
124         zd_ika = MAX( MIN(zd_ika, zco32m_br), -c_i_bio(jn_ika,layer) ) ! calcification cannot exceed CO32- stock
125                                                                        ! dissolution cannot exceed ikaite stock
126         zd_ika = zd_ika * e_i_bio(layer)                               ! bulk calcification
127
128         ! Update ikaite, dic, alk and calcium bulk concentrations
129         cbu_i_bio(jn_ika,layer) = cbu_i_bio(jn_ika,layer) + zd_ika
130         cbu_i_bio(jn_dic,layer) = cbu_i_bio(jn_dic,layer) - zd_ika
131         cbu_i_bio(jn_alk,layer) = cbu_i_bio(jn_alk,layer) - 2.*zd_ika
132         cbu_i_bio(jn_cal,layer) = cbu_i_bio(jn_cal,layer) - zd_ika
133
134         ! Retrieve brine concentrations
135         c_i_bio(jn_ika,layer) = cbu_i_bio(jn_ika,layer)/e_i_bio(layer)
136         c_i_bio(jn_dic,layer) = cbu_i_bio(jn_dic,layer)/e_i_bio(layer)
137         c_i_bio(jn_alk,layer) = cbu_i_bio(jn_alk,layer)/e_i_bio(layer)
138         c_i_bio(jn_cal,layer) = cbu_i_bio(jn_cal,layer)/e_i_bio(layer)
139
140         ! diagnostic rate of formation and saturation state
141         ika_rate(layer)  = zd_ika / ddtb
142         ika_omega(layer) = zomega
143
144         WRITE(numout,*) ' layer, zd_ika : ', layer, zd_ika
145
146      END DO
147     
148!
149!------------------------------------------------------------------------------!
150! 5) End of the routine
151!------------------------------------------------------------------------------!
152!
153
154      WRITE(numout,*) 
155      WRITE(numout,*) ' ** After CaCO3 precipitation : '
156      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
157      WRITE(numout,*)
158
159       DO jn = 1, ntra_bio
160
161         IF ( ( jn .EQ. jn_dic ) .OR.
162     &        ( jn .EQ. jn_alk ) .OR.
163     &        ( jn .EQ. jn_cal ) .OR.
164     &        ( jn .EQ. jn_ika ) ) THEN
165            WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
166            WRITE(numout,*) ' c_i_bio  : ', ( c_i_bio(jn,layer),    ! concentration of tracers in brines (mmol m-3)
167     &                                        layer = 1, nlay_bio )
168            WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer),  ! concentration of tracers in bulk ice
169     &                                        layer = 1, nlay_bio )
170         ENDIF
171
172      END DO
173
174      WRITE(numout,*)
175      WRITE(numout,*) ' End of ice_ikaite '
176      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
177     
178
179!==============================================================================|
180! End of ice_ikaite
181
182      WRITE(numout,*)
183      END
Note: See TracBrowser for help on using the repository browser.