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_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

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