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

source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 10047

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

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