New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
carb_chem.F90 in branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 9309

Last change on this file since 9309 was 9309, checked in by jpalmier, 6 years ago

JPALM -- make clearer in the code and output which atm xCO2 source is used
File size: 12.1 KB
Line 
1MODULE carb_chem_mod
2   !!======================================================================
3   !!                         ***  MODULE carb_chem_mod  ***
4   !! Calculate the carbon chemistry for the whole ocean
5   !!======================================================================
6   !! History :
7   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
8   !!----------------------------------------------------------------------
9#if defined key_roam
10   !!----------------------------------------------------------------------
11   !!                                                   key_roam?
12   !!----------------------------------------------------------------------
13
14   IMPLICIT NONE
15   PRIVATE
16     
17   PUBLIC   carb_chem        ! Called in trcbio_medusa.F90
18
19   !!----------------------------------------------------------------------
20   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
21   !! $Id$
22   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
23   !!----------------------------------------------------------------------
24
25CONTAINS
26
27   SUBROUTINE carb_chem( kt )
28      !!---------------------------------------------------------------------
29      !!                     ***  ROUTINE carb_chem  ***
30      !! This called from TRC_BIO_MEDUSA and
31      !!  - ...
32      !!----------------------------------------------------------------------
33      USE bio_medusa_mod,    ONLY: f_co2flux, f_co3, f_dcf,               &
34                                   f_h2co3, f_hco3, f_henry,              &
35                                   f_kw660, f_omarg, f_omcal,             &
36                                   f_pco2atm, f_pco2w, f_ph, f_pp0,       &
37                                   f_TALK, f_TDIC, f_xco2a,               &
38# if defined key_mocsy
39                                   zpho,                                  &
40# endif
41                                   zalk, zdic, zsal, zsil, ztmp 
42      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_0, gdepw_n,    &
43                                   gphit, mbathy, tmask
44      USE in_out_manager,    ONLY: lwp, numout
45      USE oce,               ONLY: PCO2a_in_cpl, tsb, tsn
46      USE par_kind,          ONLY: wp
47      USE par_medusa,        ONLY: jpalk, jpdic, jpdin, jpsil 
48      USE par_oce,           ONLY: jp_sal, jp_tem, jpi, jpim1, jpj,       &
49                                   jpjm1, jpk
50      USE sbc_oce,           ONLY: lk_oasis
51      USE sms_medusa,        ONLY: f2_ccd_arg, f2_ccd_cal, f3_co3,        &
52                                   f3_h2co3, f3_hco3, f3_omarg,           &
53                                   f3_omcal, f3_pH
54      USE trc,               ONLY: trn
55
56# if defined key_mocsy
57      USE mocsy_wrapper,     ONLY: mocsy_interface
58# else
59      USE trcco2_medusa,     ONLY: trc_co2_medusa
60# endif
61
62   !!* Substitution
63#  include "domzgr_substitute.h90"
64 
65      !! time (integer timestep)
66      INTEGER, INTENT( in ) ::    kt
67
68      !! Flags to help with calculating the position of the CCD
69      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
70
71      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
72      REAL(wp) :: f_rhosw
73      !! Output arguments from mocsy_interface, which aren't used
74      REAL(wp) :: f_fco2w_dum, f_BetaD_dum, f_opres_dum
75      REAL(wp) :: f_insitut_dum, f_fco2atm_dum
76      REAL(wp) :: f_schmidtco2_dum, f_kwco2_dum, f_K0_dum
77      REAL(wp) :: f_co2starair_dum, f_dpco2_dum
78      !! temporary variables
79      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
80
81      INTEGER :: iters
82      !! Loop variables
83      INTEGER :: ji, jj, jk
84
85      !!----------------------------------------------------------------------
86      !! Calculate the carbonate chemistry for the whole ocean on the first
87      !! simulation timestep and every month subsequently; the resulting 3D
88      !! field of omega calcite is used to determine the depth of the CCD
89      !!----------------------------------------------------------------------
90      !!
91!      IF(lwp) WRITE(numout,*)                                               &
92!              ' MEDUSA calculating all carbonate chemistry at kt =', kt
93!         CALL flush(numout)
94      !! blank flags
95      i2_omcal(:,:) = 0
96      i2_omarg(:,:) = 0
97
98      !! Loop over levels
99      DO jk = 1,jpk
100
101         DO jj = 2,jpjm1
102            DO ji = 2,jpim1
103               !! OPEN wet point IF..THEN loop
104               IF (tmask(ji,jj,jk).eq.1) THEN
105                  !! Do carbonate chemistry
106                  !!
107                  !! Set up required state variables.
108                  !! dissolved inorganic carbon
109                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
110                  !! alkalinity
111                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
112                  !! temperature
113                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
114                  !! salinity
115                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
116#  if defined key_mocsy
117                  !! silicic acid
118                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
119                  !! phosphate via DIN and Redfield
120                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
121#  endif
122        !!
123        !! AXY (28/02/14): check input fields
124        if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
125                     IF (lwp) WRITE(numout,*) ' carb_chem: T WARNING 3D, ',  &
126                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
127                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
128                IF(lwp) WRITE(numout,*)                                 &
129                        ' carb_chem: T SWITCHING 3D, ',                      &
130         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
131                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
132                  endif
133        if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
134                     IF (lwp) WRITE(numout,*) ' carb_chem: S WARNING 3D, ',  &
135                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
136                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
137                  endif
138               ENDIF
139            ENDDO
140         ENDDO
141
142         !! Blank input variables not used at this stage (they
143         !! relate to air-sea flux)
144         f_kw660(:,:) = 1.0
145         f_pp0(:,:)   = 1.0
146
147         DO jj = 2,jpjm1
148            DO ji = 2,jpim1
149               IF (tmask(ji,jj,jk).eq.1) THEN
150                  !! calculate carbonate chemistry at grid cell midpoint
151#  if defined key_mocsy
152                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
153                  !!                 chemistry package
154                  CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
155                                       zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
156                                       f_pp0(ji,jj),fsdept(ji,jj,jk),       &
157                                       gphit(ji,jj),f_kw660(ji,jj),         &
158                                       f_xco2a(ji,jj),1,f_ph(ji,jj),        &
159                                       f_pco2w(ji,jj),f_fco2w_dum,          &
160                                       f_h2co3(ji,jj),f_hco3(ji,jj),        &
161                                       f_co3(ji,jj),f_omarg(ji,jj),         &
162                                       f_omcal(ji,jj),f_BetaD_dum,          &
163                                       f_rhosw,f_opres_dum,                 &
164                                       f_insitut_dum,f_pco2atm(ji,jj),      &
165                                       f_fco2atm_dum,f_schmidtco2_dum,      &
166                                       f_kwco2_dum,f_K0_dum,                &
167                                       f_co2starair_dum,f_co2flux(ji,jj),   & 
168                                       f_dpco2_dum)
169                  !!
170                  !! mmol / m3 -> umol / kg
171                  f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw) * 1000.
172                  !! meq / m3 -> ueq / kg
173                  f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw) * 1000.
174                  f_dcf(ji,jj)  = f_rhosw
175#  else
176                  !! AXY (22/06/15): use old PML carbonate chemistry
177                  !! package (the MEDUSA-2 default)
178                  CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),  &
179                                      zalk(ji,jj),fsdept(ji,jj,jk),         &
180                                      f_kw660(ji,jj),f_xco2a(ji,jj),        &
181                                      f_ph(ji,jj),                          &
182                                      f_pco2w(ji,jj),f_h2co3(ji,jj),        &
183                                      f_hco3(ji,jj),f_co3(ji,jj),           &
184                                      f_omcal(ji,jj),f_omarg(ji,jj),        &
185                                      f_co2flux(ji,jj),f_TDIC(ji,jj),       &
186                                      f_TALK(ji,jj),f_dcf(ji,jj),           &
187                                      f_henry(ji,jj),iters)
188                  !!
189                  !! AXY (28/02/14): check output fields
190                  IF (iters .eq. 25) THEN
191                     IF(lwp) WRITE(numout,*)                                &
192                        ' carb_chem: 3D ITERS WARNING, ', iters, ' AT (',   &
193                        ji, ', ', jj, ', ', jk, ') AT ', kt
194                  ENDIF
195#  endif
196                  !!
197                  !! store 3D outputs
198                  f3_pH(ji,jj,jk)    = f_ph(ji,jj)
199                  f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj)
200                  f3_hco3(ji,jj,jk)  = f_hco3(ji,jj)
201                  f3_co3(ji,jj,jk)   = f_co3(ji,jj)
202                  f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
203                  f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
204                  !!
205                  !! CCD calculation: calcite
206                  if (i2_omcal(ji,jj) == 0 .and. f_omcal(ji,jj) < 1.0) then
207                     if (jk .eq. 1) then
208                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk)
209                     else
210                        fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
211                        fq1 = f3_omcal(ji,jj,jk-1) - 1.0
212                        fq2 = fq1 / (fq0 + tiny(fq0))
213                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
214                        fq4 = fq2 * fq3
215                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
216                     endif
217                     i2_omcal(ji,jj)   = 1
218                  endif
219                  if ( i2_omcal(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
220                     !! reached seafloor and still no dissolution; set
221                     !! to seafloor (W-point)
222                     f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
223                     i2_omcal(ji,jj)   = 1
224                  endif
225                  !!
226                  !! CCD calculation: aragonite
227                  if (i2_omarg(ji,jj) == 0 .and. f_omarg(ji,jj) < 1.0) then
228                     if (jk .eq. 1) then
229                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk)
230                     else
231                        fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
232                        fq1 = f3_omarg(ji,jj,jk-1) - 1.0
233                        fq2 = fq1 / (fq0 + tiny(fq0))
234                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
235                        fq4 = fq2 * fq3
236                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
237                     endif
238                     i2_omarg(ji,jj)   = 1
239                  endif
240                  if ( i2_omarg(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
241                     !! reached seafloor and still no dissolution; set
242                     !! to seafloor (W-point)
243                     f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
244                     i2_omarg(ji,jj)   = 1
245                  endif
246               ENDIF
247            ENDDO
248         ENDDO
249      ENDDO
250
251   END SUBROUTINE carb_chem
252
253#else
254   !!======================================================================
255   !!  Dummy module :                                   None key_roam?
256   !!======================================================================
257CONTAINS
258   SUBROUTINE carb_chem( )                    ! Empty routine
259      WRITE(*,*) 'carb_chem: You should not have seen this print! error?'
260   END SUBROUTINE carb_chem
261#endif 
262
263   !!======================================================================
264END MODULE carb_chem_mod
Note: See TracBrowser for help on using the repository browser.