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

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

Tidying up of headers for several files
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   !!----------------------------------------------------------------------
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# if defined key_mocsy
42                                   zpho,                                  &
43# endif
44                                   zalk, zdic, zsal, zsil, ztmp 
45      USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_0, gdepw_n,    &
46                                   gphit, mbathy, tmask
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      !! temporary variables
75      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4
76
77      INTEGER :: ji, jj, jk
78
79      !!----------------------------------------------------------------------
80      !! Calculate the carbonate chemistry for the whole ocean on the first
81      !! simulation timestep and every month subsequently; the resulting 3D
82      !! field of omega calcite is used to determine the depth of the CCD
83      !!----------------------------------------------------------------------
84      !!
85!      IF(lwp) WRITE(numout,*)                                               &
86!              ' MEDUSA calculating all carbonate chemistry at kt =', kt
87!         CALL flush(numout)
88      !! blank flags
89      i2_omcal(:,:) = 0
90      i2_omarg(:,:) = 0
91
92      !! Loop over levels
93      DO jk = 1,jpk
94
95         DO jj = 2,jpjm1
96            DO ji = 2,jpim1
97               !! OPEN wet point IF..THEN loop
98               IF (tmask(ji,jj,jk).eq.1) THEN
99                  IF (lk_oasis) THEN
100                     !! use 2D atm xCO2 from atm coupling
101                     f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
102                  ENDIF
103                  !! Do carbonate chemistry
104                  !!
105                  !! Set up required state variables.
106                  !! dissolved inorganic carbon
107                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic))
108                  !! alkalinity
109                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk))
110                  !! temperature
111                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem)
112                  !! salinity
113                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal)
114#  if defined key_mocsy
115                  !! silicic acid
116                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil))
117                  !! phosphate via DIN and Redfield
118                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
119#  endif
120        !!
121        !! AXY (28/02/14): check input fields
122        if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then
123                     IF (lwp) WRITE(numout,*) ' carb_chem: T WARNING 3D, ',  &
124                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          &
125                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
126                IF(lwp) WRITE(numout,*)                                 &
127                        ' carb_chem: T SWITCHING 3D, ',                      &
128         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem)
129                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem)     !! temperature
130                  endif
131        if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then
132                     IF (lwp) WRITE(numout,*) ' carb_chem: S WARNING 3D, ',  &
133                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          &
134                        ' at (', ji, ',', jj, ',', jk, ') at time', kt
135                  endif
136               ENDIF
137            ENDDO
138         ENDDO
139
140         !! Blank input variables not used at this stage (they
141         !! relate to air-sea flux)
142         f_kw660(:,:) = 1.0
143         f_pp0(:,:)   = 1.0
144
145         DO jj = 2,jpjm1
146            DO ji = 2,jpim1
147               IF (tmask(ji,jj,jk).eq.1) THEN
148                  !! calculate carbonate chemistry at grid cell midpoint
149#  if defined key_mocsy
150                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
151                  !!                 chemistry package
152                  CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
153                                       zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
154                                       f_pp0(ji,jj),fsdept(ji,jj,jk),       &
155                                       gphit(ji,jj),f_kw660(ji,jj),         &
156                                       f_xco2a(ji,jj),1,f_ph(ji,jj),        &
157                                       f_pco2w(ji,jj),f_fco2w(ji,jj),       &
158                                       f_h2co3(ji,jj),f_hco3(ji,jj),        &
159                                       f_co3(ji,jj),f_omarg(ji,jj),         &
160                                       f_omcal(ji,jj),f_BetaD(ji,jj),       &
161                                       f_rhosw(ji,jj),f_opres(ji,jj),       &
162                                       f_insitut(ji,jj),f_pco2atm(ji,jj),   &
163                                       f_fco2atm(ji,jj),f_schmidtco2(ji,jj),&
164                                       f_kwco2(ji,jj),f_K0(ji,jj),          &
165                                       f_co2starair(ji,jj),f_co2flux(ji,jj),& 
166                                       f_dpco2(ji,jj))
167                  !!
168                  !! mmol / m3 -> umol / kg
169                  f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000.
170                  !! meq / m3 -> ueq / kg
171                  f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
172                  f_dcf(ji,jj)  = f_rhosw(ji,jj)
173#  else
174                  !! AXY (22/06/15): use old PML carbonate chemistry
175                  !! package (the MEDUSA-2 default)
176                  CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),  &
177                                      zalk(ji,jj),fsdept(ji,jj,jk),         &
178                                      f_kw660(ji,jj),f_xco2a(ji,jj),        &
179                                      f_ph(ji,jj),                          &
180                                      f_pco2w(ji,jj),f_h2co3(ji,jj),        &
181                                      f_hco3(ji,jj),f_co3(ji,jj),           &
182                                      f_omcal(ji,jj),f_omarg(ji,jj),        &
183                                      f_co2flux(ji,jj),f_TDIC(ji,jj),       &
184                                      f_TALK(ji,jj),f_dcf(ji,jj),           &
185                                      f_henry(ji,jj),iters(ji,jj))
186                  !!
187                  !! AXY (28/02/14): check output fields
188                  IF (iters(ji,jj) .eq. 25) THEN
189                     IF(lwp) WRITE(numout,*)                                &
190                        ' carb_chem: 3D ITERS WARNING, ',                   &
191                        iters(ji,jj), ' AT (', ji, ', ', jj, ', ',          &
192                        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.