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

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

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