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

source: branches/UKMO/dev_r5518_GO6_fix_key_comp/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 9991

Last change on this file since 9991 was 9991, checked in by frrh, 6 years ago

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