source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 8441

Last change on this file since 8441 was 8441, checked in by frrh, 3 years ago

Commit changes relating to Met Office GMED ticket 339 for the modularisation of
of trcbio_medusa.F90.

Branch http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
from revisions 8394 to 8423 inclusive refer.

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   !!----------------------------------------------------------------------
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                  IF (lk_oasis) THEN
106                     !! use 2D atm xCO2 from atm coupling
107                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
108                  ENDIF
109                  !! Do carbonate chemistry
110                  !!
111                  !! Set up required state variables.
112                  !! dissolved inorganic carbon
113                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
114                  !! alkalinity
115                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
116                  !! temperature
117                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
118                  !! salinity
119                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
120#  if defined key_mocsy
121                  !! silicic acid
122                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
123                  !! phosphate via DIN and Redfield
124                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
125#  endif
126        !!
127        !! AXY (28/02/14): check input fields
128        if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
129                     IF (lwp) WRITE(numout,*) ' carb_chem: T WARNING 3D, ',  &
130                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
131                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
132                IF(lwp) WRITE(numout,*)                                 &
133                        ' carb_chem: T SWITCHING 3D, ',                      &
134         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
135                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
136                  endif
137        if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
138                     IF (lwp) WRITE(numout,*) ' carb_chem: S WARNING 3D, ',  &
139                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
140                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
141                  endif
142               ENDIF
143            ENDDO
144         ENDDO
145
146         !! Blank input variables not used at this stage (they
147         !! relate to air-sea flux)
148         f_kw660(:,:) = 1.0
149         f_pp0(:,:)   = 1.0
150
151         DO jj = 2,jpjm1
152            DO ji = 2,jpim1
153               IF (tmask(ji,jj,jk).eq.1) THEN
154                  !! calculate carbonate chemistry at grid cell midpoint
155#  if defined key_mocsy
156                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
157                  !!                 chemistry package
158                  CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
159                                       zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
160                                       f_pp0(ji,jj),fsdept(ji,jj,jk),       &
161                                       gphit(ji,jj),f_kw660(ji,jj),         &
162                                       f_xco2a(ji,jj),1,f_ph(ji,jj),        &
163                                       f_pco2w(ji,jj),f_fco2w_dum,          &
164                                       f_h2co3(ji,jj),f_hco3(ji,jj),        &
165                                       f_co3(ji,jj),f_omarg(ji,jj),         &
166                                       f_omcal(ji,jj),f_BetaD_dum,          &
167                                       f_rhosw,f_opres_dum,                 &
168                                       f_insitut_dum,f_pco2atm(ji,jj),      &
169                                       f_fco2atm_dum,f_schmidtco2_dum,      &
170                                       f_kwco2_dum,f_K0_dum,                &
171                                       f_co2starair_dum,f_co2flux(ji,jj),   & 
172                                       f_dpco2_dum)
173                  !!
174                  !! mmol / m3 -> umol / kg
175                  f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw) * 1000.
176                  !! meq / m3 -> ueq / kg
177                  f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw) * 1000.
178                  f_dcf(ji,jj)  = f_rhosw
179#  else
180                  !! AXY (22/06/15): use old PML carbonate chemistry
181                  !! package (the MEDUSA-2 default)
182                  CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),  &
183                                      zalk(ji,jj),fsdept(ji,jj,jk),         &
184                                      f_kw660(ji,jj),f_xco2a(ji,jj),        &
185                                      f_ph(ji,jj),                          &
186                                      f_pco2w(ji,jj),f_h2co3(ji,jj),        &
187                                      f_hco3(ji,jj),f_co3(ji,jj),           &
188                                      f_omcal(ji,jj),f_omarg(ji,jj),        &
189                                      f_co2flux(ji,jj),f_TDIC(ji,jj),       &
190                                      f_TALK(ji,jj),f_dcf(ji,jj),           &
191                                      f_henry(ji,jj),iters)
192                  !!
193                  !! AXY (28/02/14): check output fields
194                  IF (iters .eq. 25) THEN
195                     IF(lwp) WRITE(numout,*)                                &
196                        ' carb_chem: 3D ITERS WARNING, ', iters, ' AT (',   &
197                        ji, ', ', jj, ', ', jk, ') AT ', kt
198                  ENDIF
199#  endif
200                  !!
201                  !! store 3D outputs
202                  f3_pH(ji,jj,jk)    = f_ph(ji,jj)
203                  f3_h2co3(ji,jj,jk) = f_h2co3(ji,jj)
204                  f3_hco3(ji,jj,jk)  = f_hco3(ji,jj)
205                  f3_co3(ji,jj,jk)   = f_co3(ji,jj)
206                  f3_omcal(ji,jj,jk) = f_omcal(ji,jj)
207                  f3_omarg(ji,jj,jk) = f_omarg(ji,jj)
208                  !!
209                  !! CCD calculation: calcite
210                  if (i2_omcal(ji,jj) == 0 .and. f_omcal(ji,jj) < 1.0) then
211                     if (jk .eq. 1) then
212                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk)
213                     else
214                        fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj)
215                        fq1 = f3_omcal(ji,jj,jk-1) - 1.0
216                        fq2 = fq1 / (fq0 + tiny(fq0))
217                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
218                        fq4 = fq2 * fq3
219                        f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4
220                     endif
221                     i2_omcal(ji,jj)   = 1
222                  endif
223                  if ( i2_omcal(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
224                     !! reached seafloor and still no dissolution; set
225                     !! to seafloor (W-point)
226                     f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1)
227                     i2_omcal(ji,jj)   = 1
228                  endif
229                  !!
230                  !! CCD calculation: aragonite
231                  if (i2_omarg(ji,jj) == 0 .and. f_omarg(ji,jj) < 1.0) then
232                     if (jk .eq. 1) then
233                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk)
234                     else
235                        fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj)
236                        fq1 = f3_omarg(ji,jj,jk-1) - 1.0
237                        fq2 = fq1 / (fq0 + tiny(fq0))
238                        fq3 = fsdept(ji,jj,jk) - fsdept(ji,jj,jk-1)
239                        fq4 = fq2 * fq3
240                        f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4
241                     endif
242                     i2_omarg(ji,jj)   = 1
243                  endif
244                  if ( i2_omarg(ji,jj) == 0 .and. jk == mbathy(ji,jj) ) then
245                     !! reached seafloor and still no dissolution; set
246                     !! to seafloor (W-point)
247                     f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1)
248                     i2_omarg(ji,jj)   = 1
249                  endif
250               ENDIF
251            ENDDO
252         ENDDO
253      ENDDO
254
255   END SUBROUTINE carb_chem
256
257#else
258   !!======================================================================
259   !!  Dummy module :                                   None key_roam?
260   !!======================================================================
261CONTAINS
262   SUBROUTINE carb_chem( )                    ! Empty routine
263      WRITE(*,*) 'carb_chem: You should not have seen this print! error?'
264   END SUBROUTINE carb_chem
265#endif 
266
267   !!======================================================================
268END MODULE carb_chem_mod
Note: See TracBrowser for help on using the repository browser.