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/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 7938

Last change on this file since 7938 was 7938, checked in by marc, 7 years ago

Carbon chemistry has been pulled out of trcbio_medusa.F90 into carb_chem.F90
File size: 12.0 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: iters, f_BetaD, f_dcf, f_fco2atm,      &
34                                   f_co2flux, f_co2starair, f_co3,        &
35                                   f_fco2w, f_dpco2, f_h2co3, f_hco3,     &
36                                   f_insitut, f_K0, f_kw660, f_kwco2,     &
37                                   f_omarg, f_omcal, f_opres, f_pco2atm,  &
38                                   f_pco2w, f_ph, f_pp0, f_rhosw,         &
39                                   f_schmidtco2, f_TALK, f_TDIC, f_xco2a, &
40                                   zalk, zdic, zpho, zsal, zsil, ztmp 
41      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_n, gphit,      &
42                                   mbathy, tmask
43      USE in_out_manager,    ONLY: lwp, numout
44      USE mocsy_wrapper,     ONLY: mocsy_interface
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      USE trcco2_medusa,     ONLY: trc_co2_medusa
56
57   !!* Substitution
58#  include "domzgr_substitute.h90"
59 
60      !! time (integer timestep)
61      INTEGER, INTENT( in ) ::    kt
62
63      !! Flags to help with calculating the position of the CCD
64      INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg
65
66      !! temporary variables
67      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
68
69      INTEGER :: ji, jj, jk
70
71      !!----------------------------------------------------------------------
72      !! Calculate the carbonate chemistry for the whole ocean on the first
73      !! simulation timestep and every month subsequently; the resulting 3D
74      !! field of omega calcite is used to determine the depth of the CCD
75      !!----------------------------------------------------------------------
76      !!
77!      IF(lwp) WRITE(numout,*)                                               &
78!              ' MEDUSA calculating all carbonate chemistry at kt =', kt
79!         CALL flush(numout)
80      !! blank flags
81      i2_omcal(:,:) = 0
82      i2_omarg(:,:) = 0
83
84      !! loop over 3D space
85      DO jk = 1,jpk
86         DO jj = 2,jpjm1
87            DO ji = 2,jpim1
88               !! OPEN wet point IF..THEN loop
89               if (tmask(ji,jj,jk).eq.1) then
90                  IF (lk_oasis) THEN
91                     !! use 2D atm xCO2 from atm coupling
92                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
93                  ENDIF
94                  !! Do carbonate chemistry
95                  !!
96                  !! Set up required state variables.
97                  !! dissolved inorganic carbon
98                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
99                  !! alkalinity
100                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
101                  !! temperature
102                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
103                  !! salinity
104                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
105#  if defined key_mocsy
106                  !! silicic acid
107                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
108                  !! phosphate via DIN and Redfield
109                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
110#  endif
111        !!
112        !! AXY (28/02/14): check input fields
113        if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
114                     IF (lwp) WRITE(numout,*) ' carb_chem: T WARNING 3D, ',  &
115                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
116                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
117                IF(lwp) WRITE(numout,*)                                 &
118                        ' carb_chem: T SWITCHING 3D, ',                 &
119         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
120                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
121                  endif
122        if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
123                     IF (lwp) WRITE(numout,*) ' carb_chem: S WARNING 3D, ',  &
124                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
125                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
126                  endif
127! MAYBE BUT A BREAK IN HERE - marc
128                  !!
129                  !! Blank input variables not used at this stage (they
130                  !! relate to air-sea flux)
131                  f_kw660(ji,jj) = 1.0
132                  f_pp0(ji,jj)   = 1.0
133                  !!
134                  !! calculate carbonate chemistry at grid cell midpoint
135#  if defined key_mocsy
136                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
137                  !!                 chemistry package
138                  CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
139                                       zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
140                                       f_pp0(ji,jj),fsdept(ji,jj,jk),       &
141                                       gphit(ji,jj),f_kw660(ji,jj),         &
142                                       f_xco2a(ji,jj),1,f_ph(ji,jj),        &
143                                       f_pco2w(ji,jj),f_fco2w(ji,jj),       &
144                                       f_h2co3(ji,jj),f_hco3(ji,jj),        &
145                                       f_co3(ji,jj),f_omarg(ji,jj),         &
146                                       f_omcal(ji,jj),f_BetaD(ji,jj),       &
147                                       f_rhosw(ji,jj),f_opres(ji,jj),       &
148                                       f_insitut(ji,jj),f_pco2atm(ji,jj),   &
149                                       f_fco2atm(ji,jj),f_schmidtco2(ji,jj),&
150                                       f_kwco2(ji,jj),f_K0(ji,jj),          &
151                                       f_co2starair(ji,jj),f_co2flux(ji,jj),& 
152                                       f_dpco2(ji,jj))
153                  !!
154                  !! mmol / m3 -> umol / kg
155                  f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000.
156                  !! meq / m3 -> ueq / kg
157                  f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
158                  f_dcf(ji,jj)  = f_rhosw(ji,jj)
159#  else
160                  !! AXY (22/06/15): use old PML carbonate chemistry
161                  !! package (the MEDUSA-2 default)
162                  CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),  &
163                                      zalk(ji,jj),fsdept(ji,jj,jk),         &
164                                      f_kw660(ji,jj),f_xco2a(ji,jj),        &
165                                      f_ph(ji,jj),                          &
166                                      f_pco2w(ji,jj),f_h2co3(ji,jj),        &
167                                      f_hco3(ji,jj),f_co3(ji,jj),           &
168                                      f_omcal(ji,jj),f_omarg(ji,jj),        &
169                                      f_co2flux(ji,jj),f_TDIC(ji,jj),       &
170                                      f_TALK(ji,jj), f_dcf(ji,jj),          &
171                                      f_henry(ji,jj), iters(ji,jj))
172                  !!
173                  !! AXY (28/02/14): check output fields
174                  IF (iters(ji,jj) .eq. 25) THEN
175                     IF(lwp) WRITE(numout,*)                                &
176                        ' carb_chem: 3D ITERS WARNING, ',                   &
177                        iters(ji,jj), ' AT (', ji, ', ', jj, ', ',          &
178                        jk, ') AT ', kt
179                  ENDIF
180#  endif
181                  !!
182                  !! store 3D outputs
183                  f3_pH(ji,jj,jk)    = f_ph(ji,jj)
184                  f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj)
185                  f3_hco3(ji,jj,jk)  = f_hco3(ji,jj)
186                  f3_co3(ji,jj,jk)   = f_co3(ji,jj)
187                  f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
188                  f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
189                  !!
190                  !! CCD calculation: calcite
191                  if (i2_omcal(ji,jj) == 0 .and. f_omcal(ji,jj) < 1.0) then
192                     if (jk .eq. 1) then
193                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk)
194                     else
195                        fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
196                        fq1 = f3_omcal(ji,jj,jk-1) - 1.0
197                        fq2 = fq1 / (fq0 + tiny(fq0))
198                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
199                        fq4 = fq2 * fq3
200                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
201                     endif
202                     i2_omcal(ji,jj)   = 1
203                  endif
204                  if ( i2_omcal(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
205                     !! reached seafloor and still no dissolution; set
206                     !! to seafloor (W-point)
207                     f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
208                     i2_omcal(ji,jj)   = 1
209                  endif
210                  !!
211                  !! CCD calculation: aragonite
212                  if (i2_omarg(ji,jj) == 0 .and. f_omarg(ji,jj) < 1.0) then
213                     if (jk .eq. 1) then
214                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk)
215                     else
216                        fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
217                        fq1 = f3_omarg(ji,jj,jk-1) - 1.0
218                        fq2 = fq1 / (fq0 + tiny(fq0))
219                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
220                        fq4 = fq2 * fq3
221                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
222                     endif
223                     i2_omarg(ji,jj)   = 1
224                  endif
225                  if ( i2_omarg(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
226                     !! reached seafloor and still no dissolution; set
227                     !! to seafloor (W-point)
228                     f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
229                     i2_omarg(ji,jj)   = 1
230                  endif
231               endif
232            ENDDO
233         ENDDO
234      ENDDO
235
236   END SUBROUTINE carb_chem
237
238#else
239   !!======================================================================
240   !!  Dummy module :                                   None key_roam?
241   !!======================================================================
242CONTAINS
243   SUBROUTINE carb_chem( )                    ! Empty routine
244      WRITE(*,*) 'carb_chem: You should not have seen this print! error?'
245   END SUBROUTINE carb_chem
246#endif 
247
248   !!======================================================================
249END MODULE carb_chem_mod
Note: See TracBrowser for help on using the repository browser.