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

Last change on this file since 10149 was 10020, checked in by marc, 2 years ago

GMED ticket 406. CPP key fixes.

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