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 @ 7975

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

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