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.
air_sea.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/air_sea.F90 @ 10196

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

add DMS flux --

File size: 40.0 KB
Line 
1MODULE air_sea_mod
2   !!======================================================================
3   !!                         ***  MODULE air_sea_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   !!   -   ! 2017-08 (A. Yool)            Add air-sea flux kill switch
9   !!   -   ! 2018-08 (A. Yool)            add OMIP preindustrial DIC
10   !!   -   ! 2018-10 (A. Yool)            Add air-sea DMS flux
11   !!----------------------------------------------------------------------
12#if defined key_medusa
13   !!----------------------------------------------------------------------
14   !!                                                   MEDUSA bio-model
15   !!----------------------------------------------------------------------
16
17   IMPLICIT NONE
18   PRIVATE
19     
20   PUBLIC   air_sea        ! Called in trcbio_medusa.F90
21
22   !!----------------------------------------------------------------------
23   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
24   !! $Id$
25   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
26   !!----------------------------------------------------------------------
27
28CONTAINS
29
30   SUBROUTINE air_sea( kt )
31      !!---------------------------------------------------------------------
32      !!                     ***  ROUTINE air_sea  ***
33      !! This called from TRC_BIO_MEDUSA and
34      !!  - calculate air-sea gas exchange
35      !!  - river inputs
36      !!----------------------------------------------------------------------
37      USE bio_medusa_mod,    ONLY: f_riv_alk, f_riv_c, f_riv_n,           &
38                                   f_riv_si, f_runoff,                    & 
39                                   fgco2, zphn, zphd,                     &
40# if defined key_roam
41                                   dms_andr, dms_andr2d, dms_aran,        &
42                                   dms_aran2d, dms_hall, dms_hall2d,      &
43                                   dms_simo, dms_simo2d, dms_surf,        &
44                                   dms_surf2d, dms_andm, dms_andm2d,      &
45                                   dms_nlim, dms_wtkn,                    &
46               dms_flux, dms_flux2d,                  &
47                                   f_co2flux, f_co2flux2d,                &
48                                   f_co2starair_2d, f_co3,                &
49                                   f_dcf, f_fco2a_2d, f_fco2w_2d,         &
50                                   f_h2co3, f_hco3, f_henry,              &
51                                   f_kw660, f_kw6602d,                    &
52                                   f_o2flux, f_o2flux2d, f_o2sat,         &
53                                   f_o2sat2d, f_ocndpco2_2d,              &
54                                   f_ocnk0_2d, f_ocnkwco2_2d,             &
55                                   f_ocnrhosw_2d, f_ocnschco2_2d,         &
56                                   f_omarg, f_omcal,                      &
57                                   f_pco2a2d, f_pco2atm, f_pco2w,         &
58                                   f_pco2w2d, f_ph, f_pp0, f_pp02d,       &
59                                   f_TALK, f_TALK2d, f_TDIC, f_TDIC2d,    &
60                                   f_xco2a, f_xco2a_2d,                   &
61                                   zalk, zdic, zoxy, zsal, ztmp,          &
62#  if defined key_omip_dic
63                                   pi_fgco2,                                 &
64                                   f_pi_co2flux, f_pi_co2flux_2d,            &
65                                   f_pi_co2starair_2d, f_pi_co3,             &
66                                   f_pi_fco2a_2d, f_pi_fco2w_2d,             &
67                                   f_pi_h2co3, f_pi_hco3,                    &
68                                   f_pi_ocndpco2_2d,                         &
69                                   f_pi_omarg, f_pi_omcal,                   &
70                                   f_pi_pco2a_2d, f_pi_pco2atm, f_pi_pco2w,  &
71                                   f_pi_pco2w_2d, f_pi_ph,                   &
72                                   f_pi_TDIC, f_pi_TDIC_2d,                  &
73                                   f_pi_xco2a, f_pi_xco2a_2d,                &
74               f_pi_ph_2d, f_pi_omcal_2d, f_pi_omarg_2d, &
75               f_pi_h2co3_2d, f_pi_hco3_2d, f_pi_co3_2d, &
76                                   zomd,                                     &
77#  endif   !! key_omip_dic                               
78# endif  !! key_roam
79# if defined key_mocsy
80                                   zpho,                                  &
81# endif
82                                   zchd, zchn, zdin, zsil
83      USE dom_oce,           ONLY: e3t_0, gphit, tmask, mig, mjg
84# if defined key_vvl
85      USE dom_oce,           ONLY: e3t_n
86# endif
87      USE iom,               ONLY: lk_iomput
88      USE in_out_manager,    ONLY: lwp, numout
89      USE par_kind,          ONLY: wp
90      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
91      USE sbc_oce,           ONLY: fr_i, qsr, wndm
92      USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          &
93                                   jriver_alk, jriver_c,                  &
94                                   jriver_n, jriver_si,                   &
95                                   riv_alk, riv_c, riv_n, riv_si,         &
96                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    &
97                                   zn_dms_mld, zn_dms_qsr,                &
98                                   xnln, xnld 
99      USE trc,               ONLY: med_diag
100      USE zdfmxl,            ONLY: hmld
101
102# if defined key_roam
103      USE gastransfer,       ONLY: gas_transfer
104#  if defined key_mocsy
105      USE mocsy_wrapper,     ONLY: mocsy_interface
106#  else
107      USE trcco2_medusa,     ONLY: trc_co2_medusa
108#  endif
109      USE trcdms_medusa,     ONLY: trc_dms_medusa, dms_flux_ocn
110      USE trcoxy_medusa,     ONLY: trc_oxy_medusa
111# endif
112      USE lib_mpp,           ONLY: ctl_stop
113      USE trcstat,           ONLY: trc_rst_dia_stat 
114
115   !!* Substitution
116#  include "domzgr_substitute.h90"
117
118      !! time (integer timestep)
119      INTEGER, INTENT( in ) :: kt
120
121      !! Loop variables
122      INTEGER :: ji, jj
123
124# if defined key_roam
125      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
126      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
127      REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d
128      REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day
129
130      INTEGER :: iters
131
132      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
133      REAL(wp), DIMENSION(jpi,jpj) :: f_fco2w, f_rhosw
134      REAL(wp), DIMENSION(jpi,jpj) :: f_fco2atm
135      REAL(wp), DIMENSION(jpi,jpj) :: f_schmidtco2, f_kwco2, f_K0
136      REAL(wp), DIMENSION(jpi,jpj) :: f_co2starair, f_dpco2
137      !! Output arguments from mocsy_interface, which aren't used
138      REAL(wp) :: f_BetaD_dum, f_opres_dum
139      REAL(wp) :: f_insitut_dum
140      REAL(wp) :: f_kwo2_dum
141#  if defined key_omip_dic
142      !! AXY (06/08/18): output variables that are not different for OMIP DIC
143      REAL(wp), DIMENSION(jpi,jpj) :: f_pi_fco2w
144      REAL(wp), DIMENSION(jpi,jpj) :: f_pi_fco2atm
145      REAL(wp), DIMENSION(jpi,jpj) :: f_pi_co2starair, f_pi_dpco2
146      REAL(wp) :: f_rhosw_dum, f_schmidtco2_dum, f_kwco2_dum,f_K0_dum
147#  endif
148# endif
149
150# if defined key_roam
151      !! init
152      f_fco2w(:,:)       = 0.0
153      f_fco2atm(:,:)     = 0.0
154      f_schmidtco2(:,:)  = 0.0
155      f_kwco2(:,:)       = 0.0
156      f_co2starair(:,:)  = 0.0
157      f_dpco2(:,:)       = 0.0
158      f_rhosw(:,:)       = 0.0
159      f_K0(:,:)          = 0.0
160#  if defined key_omip_dic
161      f_pi_fco2w(:,:)       = 0.0
162      f_pi_fco2atm(:,:)     = 0.0
163      f_pi_co2starair(:,:)  = 0.0
164      f_pi_dpco2(:,:)       = 0.0     
165#  endif
166      !! air pressure (atm); ultimately this will use air
167      !! pressure at the base of the UKESM1 atmosphere
168      !!                                     
169      f_pp0(:,:)   = 1.0
170
171      !!-----------------------------------------------------------
172      !! Air-sea gas exchange
173      !!-----------------------------------------------------------
174
175#  if defined key_debug_medusa
176               IF (lwp) write (numout,*)                     & 
177               'air-sea: gas_transfer kt = ', kt
178               CALL flush(numout)
179#  endif
180      DO jj = 2,jpjm1
181         DO ji = 2,jpim1
182            !! OPEN wet point IF..THEN loop
183            IF (tmask(ji,jj,1) == 1) then
184               !!
185               !! AXY (23/06/15): as part of an effort to update the
186               !!                 carbonate chemistry in MEDUSA, the gas
187               !!                 transfer velocity used in the carbon
188               !!                 and oxygen cycles has been harmonised
189               !!                 and is calculated by the same function
190               !!                 here; this harmonisation includes
191               !!                 changes to the PML carbonate chemistry
192               !!                 scheme so that it too makes use of the
193               !!                 same gas transfer velocity; the
194               !!                 preferred parameterisation of this is
195               !!                 Wanninkhof (2014), option 7
196               !!
197               CALL gas_transfer( wndm(ji,jj), 1, 7,         &  ! inputs
198                                  f_kw660(ji,jj) )              ! outputs
199            ENDIF
200         ENDDO
201      ENDDO
202
203#  if defined key_debug_medusa
204      IF (lwp) write (numout,*)                     &
205           'air-sea: carb-chem kt = ', kt
206      CALL flush(numout)
207      !! JPALM add carb print:
208      call trc_rst_dia_stat(f_xco2a(:,:), 'f_xco2a')
209      call trc_rst_dia_stat(wndm(:,:), 'wndm')
210      call trc_rst_dia_stat(f_kw660(:,:), 'f_kw660')
211      call trc_rst_dia_stat(ztmp(:,:), 'ztmp')
212      call trc_rst_dia_stat(zsal(:,:), 'zsal')
213      call trc_rst_dia_stat(zalk(:,:), 'zalk')
214      call trc_rst_dia_stat(zdic(:,:), 'zdic')
215      call trc_rst_dia_stat(zsil(:,:), 'zsil')
216      call trc_rst_dia_stat(zpho(:,:), 'zpho')
217#  endif
218#  if defined key_axy_carbchem
219#   if defined key_mocsy
220      !!-----------------------------------------------------------
221      !! MOCSY carbonate chemistry
222      !!-----------------------------------------------------------
223      !!
224      DO jj = 2,jpjm1
225         DO ji = 2,jpim1
226            if (tmask(ji,jj,1) == 1) then
227               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
228               !!          carbonate failure in the coupled run.
229               !!          must be associated to air-sea flux or air xCO2...
230               !!          Check MOCSY inputs
231               IF ( (zsal(ji,jj) > 75.0 ).OR.(zsal(ji,jj) < 0.0 ) .OR.        &
232                    (ztmp(ji,jj) > 50.0 ).OR.(ztmp(ji,jj) < -20.0 ) .OR.      &
233                    (zalk(ji,jj) > 35.0E2 ).OR.(zalk(ji,jj) <= 0.0 ) .OR.     &
234                    (zdic(ji,jj) > 35.0E2 ).OR.(zdic(ji,jj) <= 0.0 ) .OR.     &
235                    (f_kw660(ji,jj) > 1.0E-2 ).OR.(f_kw660(ji,jj) < 0.0 ) ) THEN
236                  IF(lwp) THEN
237                     WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
238                     WRITE(numout,*) ' surface S = ',zsal(ji,jj)
239                     WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
240                     WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
241                     WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
242                     WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
243                     WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
244                     WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
245                     WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
246                     WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
247                     WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
248                     WRITE(numout,*) ' MOCSY input: ji =', mig(ji),' jj = ', mjg(jj),  &
249                          ' kt = ', kt 
250                     WRITE(numout,*) 'MEDUSA - Air-Sea INPUT: unrealistic surface Carb. Chemistry'
251                  ENDIF
252                  CALL ctl_stop( 'MEDUSA - Air-Sea INPUT: ',             &
253                                 'unrealistic surface Carb. Chemistry -- INPUTS' )
254               ENDIF
255               !!
256               !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
257               !!                 chemistry package; note that depth is set to
258               !!                 zero in this call
259               CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj),     &
260                                    zdic(ji,jj),zsil(ji,jj),zpho(ji,jj),     &
261                                    f_pp0(ji,jj),0.0,                        &
262                                    gphit(ji,jj),f_kw660(ji,jj),             &
263                                    f_xco2a(ji,jj),1,f_ph(ji,jj),            &
264                                    f_pco2w(ji,jj),f_fco2w(ji,jj),           &
265                                    f_h2co3(ji,jj),f_hco3(ji,jj),            &
266                                    f_co3(ji,jj),f_omarg(ji,jj),             &
267                                    f_omcal(ji,jj),f_BetaD_dum,              &
268                                    f_rhosw(ji,jj),f_opres_dum,              &
269                                    f_insitut_dum,f_pco2atm(ji,jj),          &
270                                    f_fco2atm(ji,jj),f_schmidtco2(ji,jj),    &
271                                    f_kwco2(ji,jj),f_K0(ji,jj),              &
272                                    f_co2starair(ji,jj),f_co2flux(ji,jj),    &
273                                    f_dpco2(ji,jj))
274               !! mmol / m3 -> umol / kg
275               f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000.
276               !! meq / m3 ->  ueq / kg
277               f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
278               f_dcf(ji,jj)  = f_rhosw(ji,jj)
279               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
280               !!          carbonate failure in the coupled run.
281               !!          must be associated to air-sea flux or air xCO2...
282               !!          Check MOCSY outputs
283               !!===================
284               !! Jpalm -- 19-02-2018 -- remove the cap - only check MOCSY inputs
285               !!       because of specific area in arabic sea where strangely
286               !!       with core 2 forcing, ALK is lower than DIC and result in
287               !!       Enormous dpco2 - even if all carb chem caract are OK.
288               !!       and this check stops the model.
289               !!       --Input checks are already more than enough to stop the
290               !!       model if carb chem goes crazy.
291               !!       we remove the mocsy output checks
292               !!===================
293               !!IF ( (f_pco2w(ji,jj) > 1.E4 ).OR.(f_pco2w(ji,jj) < 0.0 ) .OR.     &
294               !!    (f_fco2w(ji,jj) > 1.E4 ).OR.(f_fco2w(ji,jj) < 0.0 ) .OR.     &   
295               !!    (f_fco2atm(ji,jj) > 1.E4 ).OR.(f_fco2atm(ji,jj) < 0.0 ) .OR.     &
296               !!    (f_co2flux(ji,jj) > 1.E-1 ).OR.(f_co2flux(ji,jj) < -1.E-1 ) .OR.     &
297               !!    (f_dpco2(ji,jj) > 1.E4 ).OR.(f_dpco2(ji,jj) < -1.E4 ) ) THEN
298               !!  IF(lwp) THEN
299               !!      WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
300               !!      WRITE(numout,*) ' surface S = ',zsal(ji,jj)
301               !!      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
302               !!      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
303               !!      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
304               !!      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
305               !!      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
306               !!      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
307               !!      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
308               !!      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
309               !!      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
310               !!      WRITE(numout,*) ' MOCSY output: ji =', mig(ji),' jj = ', mjg(jj),  &
311               !!                        ' kt = ', kt     
312               !!      WRITE(numout,*) 'MEDUSA - Air-Sea OUTPUT: unrealistic surface Carb. Chemistry'
313               !!  ENDIF     
314               !!  CALL ctl_stop( 'MEDUSA - Air-Sea OUTPUT: ',            &
315               !!                 'unrealistic surface Carb. Chemistry -- OUTPUTS' )
316               !!ENDIF     
317#    if defined key_omip_dic
318               !! AXY (06/08/18): repeat carbonate chemistry calculations but
319               !!                 using preindustrial DIC for OMIP; note that
320               !!                 some outputs below are dummy since they are
321               !!                 the same (or should be!) as for regular DIC
322               CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj),        &
323                                    zomd(ji,jj),zsil(ji,jj),zpho(ji,jj),        &
324                                    f_pp0(ji,jj),0.0,                           &
325                                    gphit(ji,jj),f_kw660(ji,jj),                &
326                                    f_pi_xco2a(ji,jj),1,f_pi_ph(ji,jj),         &
327                                    f_pi_pco2w(ji,jj),f_pi_fco2w(ji,jj),        &
328                                    f_pi_h2co3(ji,jj),f_pi_hco3(ji,jj),         &
329                                    f_pi_co3(ji,jj),f_pi_omarg(ji,jj),          &
330                                    f_pi_omcal(ji,jj),f_BetaD_dum,              &
331                                    f_rhosw_dum,f_opres_dum,                    &
332                                    f_insitut_dum,f_pi_pco2atm(ji,jj),          &
333                                    f_pi_fco2atm(ji,jj),f_schmidtco2_dum,       &
334                                    f_kwco2_dum,f_K0_dum,                       &
335                                    f_pi_co2starair(ji,jj),f_pi_co2flux(ji,jj), &
336                                    f_pi_dpco2(ji,jj))
337               !! mmol / m3 -> umol / kg
338               f_pi_TDIC(ji,jj) = (zomd(ji,jj) / f_rhosw(ji,jj)) * 1000.
339#    endif  !! omip
340            ENDIF
341         ENDDO
342      ENDDO
343
344#    if defined key_debug_medusa
345      !! JPALM add carb print:
346      call trc_rst_dia_stat(f_pco2w(:,:), 'f_pco2w')
347      call trc_rst_dia_stat(f_fco2w(:,:), 'f_fco2w')
348      call trc_rst_dia_stat(f_fco2atm(:,:), 'f_fco2atm')
349      call trc_rst_dia_stat(f_schmidtco2(:,:), 'f_schmidtco2')
350      call trc_rst_dia_stat(f_kwco2(:,:), 'f_kwco2')
351      call trc_rst_dia_stat(f_co2starair(:,:), 'f_co2starair')
352      call trc_rst_dia_stat(f_co2flux(:,:), 'f_co2flux')
353      call trc_rst_dia_stat(f_dpco2(:,:), 'f_dpco2')
354#     if defined key_omip_dic
355      call trc_rst_dia_stat(f_pi_pco2w(:,:), 'f_pi_pco2w')
356      call trc_rst_dia_stat(f_pi_fco2w(:,:), 'f_pi_fco2w')
357      call trc_rst_dia_stat(f_pi_fco2atm(:,:), 'f_pi_fco2atm')
358      call trc_rst_dia_stat(f_pi_co2starair(:,:), 'f_pi_co2starair')
359      call trc_rst_dia_stat(f_pi_co2flux(:,:), 'f_pi_co2flux')
360      call trc_rst_dia_stat(f_pi_dpco2(:,:), 'f_pi_dpco2')     
361#     endif
362#    endif
363               
364#   else   !! mocsy
365      !!-----------------------------------------------------------
366      !! Blackford et al. (2007) carbonate chemistry
367      !!-----------------------------------------------------------
368      !!
369      DO jj = 2,jpjm1
370         DO ji = 2,jpim1
371            if (tmask(ji,jj,1) == 1) then     
372               iters = 0
373               !!
374               !! carbon dioxide (CO2); Jerry Blackford code (ostensibly
375               !! OCMIP-2, but not); this should be removed in the future
376               CALL trc_co2_medusa(ztmp(ji,jj),zsal(ji,jj),zdic(ji,jj),      &
377                                   zalk(ji,jj),0.0,                          &
378                                   f_kw660(ji,jj),f_xco2a(ji,jj),            &
379                                   f_ph(ji,jj),                              &
380                                   f_pco2w(ji,jj),f_h2co3(ji,jj),            &
381                                   f_hco3(ji,jj),f_co3(ji,jj),               &
382                                   f_omcal(ji,jj),f_omarg(ji,jj),            &
383                                   f_co2flux(ji,jj),f_TDIC(ji,jj),           &
384                                   f_TALK(ji,jj),f_dcf(ji,jj),               &
385                                   f_henry(ji,jj),iters)
386               !!
387               !! AXY (09/01/14): removed iteration and NaN checks; these have
388               !!                 been moved to trc_co2_medusa together with a
389               !!                 fudge that amends erroneous values (this is
390               !!                 intended to be a temporary fudge!); the
391               !!                 output warnings are retained here so that
392               !!                 failure position can be determined
393               if (iters .eq. 25) then
394                  IF(lwp) WRITE(numout,*) 'air-sea: ITERS WARNING, ',       &
395                       iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
396               endif
397            ENDIF
398         ENDDO
399      ENDDO
400
401#   endif !! mocsy
402#  else  !! carbchem
403      !!-----------------------------------------------------------
404      !! No carbonate chemistry
405      !!-----------------------------------------------------------
406      !!
407      DO jj = 2,jpjm1
408         DO ji = 2,jpim1
409            if (tmask(ji,jj,1) == 1) then
410               !! AXY (18/04/13): switch off carbonate chemistry
411               !!                 calculations; provide quasi-sensible
412               !!                 alternatives
413               f_ph(ji,jj)           = 8.1
414               f_pco2w(ji,jj)        = f_xco2a(ji,jj)
415               f_h2co3(ji,jj)        = 0.005 * zdic(ji,jj)
416               f_hco3(ji,jj)         = 0.865 * zdic(ji,jj)
417               f_co3(ji,jj)          = 0.130 * zdic(ji,jj)
418               f_omcal(ji,jj) = 4.
419               f_omarg(ji,jj) = 2.
420               f_co2flux(ji,jj)      = 0.
421               f_TDIC(ji,jj)         = zdic(ji,jj)
422               f_TALK(ji,jj)         = zalk(ji,jj)
423               f_dcf(ji,jj)          = 1.026
424               f_henry(ji,jj)        = 1.
425               !! AXY (23/06/15): add in some extra MOCSY diagnostics
426               f_fco2w(ji,jj)        = f_xco2a(ji,jj)
427! This doesn't seem to be used - marc 16/5/17
428!               f_BetaD(ji,jj)        = 1.
429               f_rhosw(ji,jj)        = 1.026
430! This doesn't seem to be used - marc 16/5/17
431!               f_opres(ji,jj)        = 0.
432!               f_insitut(ji,jj)      = ztmp(ji,jj)
433               f_pco2atm(ji,jj)      = f_xco2a(ji,jj)
434               f_fco2atm(ji,jj)      = f_xco2a(ji,jj)
435               f_schmidtco2(ji,jj)   = 660.
436               f_kwco2(ji,jj)        = 0.
437               f_K0(ji,jj)           = 0.
438               f_co2starair(ji,jj)   = f_xco2a(ji,jj)
439               f_dpco2(ji,jj)        = 0.
440            ENDIF
441         ENDDO
442      ENDDO
443#  endif  !! carbchem
444
445#  if defined key_axy_killco2flux
446      !! AXY (18/08/17): single kill switch on air-sea CO2 flux for budget checking
447      f_co2flux(:,:) = 0.
448#   if defined key_omip_dic
449      !! AXY (06/08/18): applies to OMIP DIC too
450      f_pi_co2flux(:,:) = 0.
451#   endif     
452#  endif
453
454      DO jj = 2,jpjm1
455         DO ji = 2,jpim1
456            if (tmask(ji,jj,1) == 1) then
457               !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide
458               !! through by layer thickness
459               f_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_co2flux(ji,jj) *       &
460                                  86400. / fse3t(ji,jj,1)
461#   if defined key_omip_dic
462               !! PI DIC - for omip runs
463               f_pi_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_pi_co2flux(ji,jj) * &
464                                      86400. / fse3t(ji,jj,1)
465#   endif
466               !! oxygen (O2); OCMIP-2 code
467               !! AXY (23/06/15): amend input list for oxygen to account
468               !!                 for common gas transfer velocity
469               CALL trc_oxy_medusa(ztmp(ji,jj),zsal(ji,jj),f_kw660(ji,jj),     &
470                                   f_pp0(ji,jj),zoxy(ji,jj),                   &
471                                   f_kwo2_dum,f_o2flux(ji,jj),                 &
472                                   f_o2sat(ji,jj))
473               !!
474               !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide
475               !! through by layer thickness
476               f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) *       &
477                                  86400. / fse3t(ji,jj,1)
478            ENDIF
479         ENDDO
480      ENDDO
481
482      !! Jpalm (08-2014)
483      !! DMS surface concentration calculation
484      !! initialy added for UKESM1 model.
485      !! using MET-OFFICE subroutine.
486      !! DMS module only needs Chl concentration and MLD
487      !! to get an aproximate value of DMS concentration.
488      !! air-sea fluxes are calculated by atmospheric chemitry model
489      !! from atm and oc-surface concentrations.
490      !!
491      !! AXY (13/03/15): this is amended to calculate all of the DMS
492      !!                 estimates examined during UKESM1 (see
493      !!                 comments in trcdms_medusa.F90)
494      !!
495      !! AXY (25/05/17): amended to additionally pass DIN limitation as well as [DIN];
496      !!                 accounts for differences in nutrient half-saturations; changes
497      !!                 also made in trc_dms_medusa; this permits an additional DMS
498      !!                 calculation while retaining the existing Anderson one
499      !!
500      IF (jdms == 1) THEN
501         DO jj = 2,jpjm1
502            DO ji = 2,jpim1
503               if (tmask(ji,jj,1) == 1) then
504                  !! calculate weighted half-saturation for DIN uptake
505                  dms_wtkn(ji,jj) = ((zphn(ji,jj) * xnln) +            &
506                                     (zphd(ji,jj) * xnld)) /           &
507                                     (zphn(ji,jj) + zphd(ji,jj))       
508                  !!
509                  !! feed in correct inputs
510                  if (jdms_input .eq. 0) then
511                     !! use instantaneous inputs
512                     dms_nlim(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + dms_wtkn(ji,jj))
513                     !!
514                     CALL trc_dms_medusa(zchn(ji,jj),zchd(ji,jj),             &
515                                         hmld(ji,jj),qsr(ji,jj),              &
516                                         zdin(ji,jj), dms_nlim(ji,jj),        &
517                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
518                                         dms_andm)
519                  else
520                     !! use diel-average inputs
521                     dms_nlim(ji,jj) = zn_dms_din(ji,jj) /                    &
522                                      (zn_dms_din(ji,jj) + dms_wtkn(ji,jj))
523                     !!
524                     CALL trc_dms_medusa(zn_dms_chn(ji,jj),zn_dms_chd(ji,jj), &
525                                         zn_dms_mld(ji,jj),zn_dms_qsr(ji,jj), &
526                                         zn_dms_din(ji,jj),dms_nlim(ji,jj),   &
527                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
528                                         dms_andm)
529                  endif
530                  !!
531                  !! assign correct output to variable passed to atmosphere
532                  if (jdms_model .eq. 1) then
533                     dms_surf = dms_andr
534                  elseif (jdms_model .eq. 2) then
535                     dms_surf = dms_simo
536                  elseif (jdms_model .eq. 3) then
537                     dms_surf = dms_aran
538                  elseif (jdms_model .eq. 4) then
539                     dms_surf = dms_hall
540                  elseif (jdms_model .eq. 5) then
541                     dms_surf = dms_andm
542                  endif
543                  !!
544                  !! calculate air-sea DMS flux using UM code
545                  CALL dms_flux_ocn(wndm(ji,jj),ztmp(ji,jj),dms_surf,1, &
546                                    dms_flux)
547                  !!
548                  !! Correct DMS flux for sea-ice, but leave units as kg S m-2 s-1
549                  dms_flux = (1. - fr_i(ji,jj)) * dms_flux
550                  !!
551                  !! 2D diag through iom_use
552                  IF( med_diag%DMS_SURF%dgsave ) THEN
553                     dms_surf2d(ji,jj) = dms_surf
554                  ENDIF
555                  IF( med_diag%DMS_ANDR%dgsave ) THEN
556                     dms_andr2d(ji,jj) = dms_andr
557                  ENDIF
558                  IF( med_diag%DMS_SIMO%dgsave ) THEN
559                     dms_simo2d(ji,jj) = dms_simo
560                  ENDIF
561                  IF( med_diag%DMS_ARAN%dgsave ) THEN
562                     dms_aran2d(ji,jj) = dms_aran
563                  ENDIF
564                  IF( med_diag%DMS_HALL%dgsave ) THEN
565                     dms_hall2d(ji,jj) = dms_hall
566                  ENDIF
567                  IF( med_diag%DMS_ANDM%dgsave ) THEN
568                     dms_andm2d(ji,jj) = dms_andm
569                  ENDIF
570                  IF( med_diag%DMS_FLUX%dgsave ) THEN
571                     dms_flux2d(ji,jj) = dms_flux
572                  ENDIF
573               ENDIF
574            ENDDO
575         ENDDO
576#   if defined key_debug_medusa
577         IF (lwp) write (numout,*)                                &
578            'air-sea: finish calculating dms kt = ',kt
579            CALL flush(numout)
580#   endif
581      ENDIF                  !! End IF (jdms == 1)
582
583      !!
584      !! store 2D outputs
585      !!
586      !! JPALM -- 17-11-16 -- put fgco2 out of diag request
587      !!       is needed for coupling; pass through restart
588      DO jj = 2,jpjm1
589         DO ji = 2,jpim1
590            if (tmask(ji,jj,1) == 1) then
591               !! IF( med_diag%FGCO2%dgsave ) THEN
592               !! convert from  mol/m2/day to kg/m2/s
593               !! mmol-C/m3/d -> kg-CO2/m2/s
594               fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) *             &
595                              CO2flux_conv
596               !! ENDIF
597               IF ( lk_iomput ) THEN
598                  IF( med_diag%ATM_PCO2%dgsave ) THEN
599                     f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
600                  ENDIF
601                  IF( med_diag%OCN_PCO2%dgsave ) THEN
602                     f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
603                  ENDIF
604                  IF( med_diag%CO2FLUX%dgsave ) THEN
605                     !! mmol/m3/d -> mmol/m2/d
606                     f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1)
607                  ENDIF
608                  IF( med_diag%TCO2%dgsave ) THEN
609                     f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
610                  ENDIF
611                  IF( med_diag%TALK%dgsave ) THEN
612                     f_TALK2d(ji,jj) = f_TALK(ji,jj)
613                  ENDIF
614                  IF( med_diag%KW660%dgsave ) THEN
615                      f_kw6602d(ji,jj) = f_kw660(ji,jj)
616                  ENDIF
617                  IF( med_diag%ATM_PP0%dgsave ) THEN
618                      f_pp02d(ji,jj) = f_pp0(ji,jj)
619                  ENDIF
620                  IF( med_diag%O2FLUX%dgsave ) THEN
621                     f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
622                  ENDIF
623                  IF( med_diag%O2SAT%dgsave ) THEN
624                     f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
625                  ENDIF
626                  !! AXY (24/11/16): add in extra MOCSY diagnostics
627                  IF( med_diag%ATM_XCO2%dgsave ) THEN
628                     f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
629                  ENDIF
630                  IF( med_diag%OCN_FCO2%dgsave ) THEN
631                     f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
632                  ENDIF
633                  IF( med_diag%ATM_FCO2%dgsave ) THEN
634                     f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
635                  ENDIF
636                  IF( med_diag%OCN_RHOSW%dgsave ) THEN
637                     f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
638                  ENDIF
639                  IF( med_diag%OCN_SCHCO2%dgsave ) THEN
640                     f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
641                  ENDIF
642                  IF( med_diag%OCN_KWCO2%dgsave ) THEN
643                     f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
644                  ENDIF
645                  IF( med_diag%OCN_K0%dgsave ) THEN
646                     f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
647                  ENDIF
648                  IF( med_diag%CO2STARAIR%dgsave ) THEN
649                     f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
650                  ENDIF
651                  IF( med_diag%OCN_DPCO2%dgsave ) THEN
652                     f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
653                  ENDIF
654#   if defined key_omip_dic
655                  !! AXY (06/08/18): diagnostics for OMIP PI DIC
656                  IF( med_diag%PI_ATM_PCO2%dgsave ) THEN
657                     f_pi_pco2a_2d(ji,jj) = f_pi_pco2atm(ji,jj)
658                  ENDIF
659                  IF( med_diag%PI_OCN_PCO2%dgsave ) THEN
660                     f_pi_pco2w_2d(ji,jj) = f_pi_pco2w(ji,jj)
661                  ENDIF
662                  IF( med_diag%PI_CO2FLUX%dgsave ) THEN
663                     !! mmol/m3/d -> mmol/m2/d
664                     f_pi_co2flux_2d(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1)
665                  ENDIF
666                  IF( med_diag%PI_TCO2%dgsave ) THEN
667                     f_PI_TDIC_2d(ji,jj) = f_PI_TDIC(ji,jj)
668                  ENDIF
669                  IF( med_diag%PI_ATM_XCO2%dgsave ) THEN
670                     f_pi_xco2a_2d(ji,jj) = f_pi_xco2a(ji,jj)
671                  ENDIF
672                  IF( med_diag%PI_OCN_FCO2%dgsave ) THEN
673                     f_pi_fco2w_2d(ji,jj) = f_pi_fco2w(ji,jj)
674                  ENDIF
675                  IF( med_diag%PI_ATM_FCO2%dgsave ) THEN
676                     f_pi_fco2a_2d(ji,jj) = f_pi_fco2atm(ji,jj)
677                  ENDIF
678                  IF( med_diag%PI_CO2STARAIR%dgsave ) THEN
679                     f_pi_co2starair_2d(ji,jj) = f_pi_co2starair(ji,jj)
680                  ENDIF
681                  IF( med_diag%PI_OCN_DPCO2%dgsave ) THEN
682                     f_pi_ocndpco2_2d(ji,jj) = f_pi_dpco2(ji,jj)
683                  ENDIF                 
684                  IF( med_diag%PI_FGCO2%dgsave ) THEN
685                     !! convert from  mol/m2/day to kg/m2/s
686                     !! mmol-C/m3/d -> kg-CO2/m2/s
687                     pi_fgco2(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1) * &
688                          CO2flux_conv
689                  ENDIF
690                  IF( med_diag%PI_OCN_PH%dgsave ) THEN
691                     f_pi_ph_2d(ji,jj) = f_pi_ph(ji,jj)
692                  ENDIF                 
693                  IF( med_diag%PI_OCNH2CO3%dgsave ) THEN
694                     f_pi_h2co3_2d(ji,jj) = f_pi_h2co3(ji,jj)
695                  ENDIF                 
696                  IF( med_diag%PI_OCN_HCO3%dgsave ) THEN
697                     f_pi_hco3_2d(ji,jj) = f_pi_hco3(ji,jj)
698                  ENDIF                 
699                  IF( med_diag%PI_OCN_CO3%dgsave ) THEN
700                     f_pi_co3_2d(ji,jj) = f_pi_co3(ji,jj)
701                  ENDIF                 
702                  IF( med_diag%PI_OM_CAL%dgsave ) THEN
703                     f_pi_omcal_2d(ji,jj) = f_pi_omcal(ji,jj)
704                  ENDIF                 
705                  IF( med_diag%PI_OM_ARG%dgsave ) THEN
706                     f_pi_omarg_2d(ji,jj) = f_pi_omarg(ji,jj)
707                  ENDIF                 
708#   endif                 
709               ENDIF
710            ENDIF
711         ENDDO
712      ENDDO
713
714      !!-----------------------------------------------------------------
715      !! River inputs
716      !!-----------------------------------------------------------------
717      DO jj = 2,jpjm1
718         DO ji = 2,jpim1
719            !! OPEN wet point IF..THEN loop
720            if (tmask(ji,jj,1) == 1) then
721               !!
722               !! runoff comes in as        kg / m2 / s
723               !! used and written out as   m3 / m2 / d (= m / d)
724               !! where                     1000 kg / m2 / d =
725               !!                             1 m3 / m2 / d = 1 m / d
726               !!
727               !! AXY (17/07/14): the compiler doesn't like this line for
728               !!                 some reason; as MEDUSA doesn't even use
729               !!                 runoff for riverine inputs, a temporary
730               !!                 solution is to switch off runoff entirely
731               !!                 here; again, this change is one of several
732               !!                 that will need revisiting once MEDUSA has
733               !!                 bedded down in UKESM1; particularly so if
734               !!                 the land scheme provides information
735               !!                 concerning nutrient fluxes
736               !!
737               !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. *   &
738               !!                   60. * 24.
739               f_runoff(ji,jj) = 0.0
740               !!
741               !! nutrients are added via rivers to the model in one of
742               !! two ways:
743               !!   1. via river concentration; i.e. the average nutrient
744               !!      concentration of a river water is described by a
745               !!      spatial file, and this is multiplied by runoff to
746               !!      give a nutrient flux
747               !!   2. via direct river flux; i.e. the average nutrient
748               !!      flux due to rivers is described by a spatial file,
749               !!      and this is simply applied as a direct nutrient
750               !!      flux (i.e. it does not relate or respond to model
751               !!      runoff) nutrient fields are derived from the
752               !!      GlobalNEWS 2 database; carbon and alkalinity are
753               !!      derived from continent-scale DIC estimates (Huang et
754               !!      al., 2012) and some Arctic river alkalinity
755               !!      estimates (Katya?)
756               !!
757               !! as of 19/07/12, riverine nutrients can now be spread
758               !! vertically across several grid cells rather than just
759               !! poured into the surface box; this block of code is still
760               !! executed, however, to set up the total amounts of
761               !! nutrient entering via rivers
762               !!
763               !! nitrogen
764               if (jriver_n .eq. 1) then
765                  !! river concentration specified; use runoff to
766                  !! calculate input
767                  f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
768               elseif (jriver_n .eq. 2) then
769                  !! river flux specified; independent of runoff
770                  f_riv_n(ji,jj) = riv_n(ji,jj)
771               endif
772               !!
773               !! silicon
774               if (jriver_si .eq. 1) then
775                  !! river concentration specified; use runoff to
776                  !! calculate input
777                  f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
778               elseif (jriver_si .eq. 2) then
779                  !! river flux specified; independent of runoff
780                  f_riv_si(ji,jj) = riv_si(ji,jj)
781               endif
782               !!
783               !! carbon
784               if (jriver_c .eq. 1) then
785                  !! river concentration specified; use runoff to
786                  !! calculate input
787                  f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
788               elseif (jriver_c .eq. 2) then
789                  !! river flux specified; independent of runoff
790                  f_riv_c(ji,jj) = riv_c(ji,jj)
791               endif
792               !!
793               !! alkalinity
794               if (jriver_alk .eq. 1) then
795                  !! river concentration specified; use runoff to
796                  !! calculate input
797                  f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
798               elseif (jriver_alk .eq. 2) then
799                  !! river flux specified; independent of runoff
800                  f_riv_alk(ji,jj) = riv_alk(ji,jj)
801               endif
802            ENDIF
803         ENDDO
804      ENDDO
805# endif  !! Roam
806     
807   END SUBROUTINE air_sea
808
809#else  !! medusa
810   !!======================================================================
811   !!  Dummy module :                                   No MEDUSA bio-model
812   !!======================================================================
813CONTAINS
814   SUBROUTINE air_sea( )                    ! Empty routine
815      WRITE(*,*) 'air_sea: You should not have seen this print! error?'
816   END SUBROUTINE air_sea
817#endif  !! medusa
818
819   !!======================================================================
820END MODULE air_sea_mod
Note: See TracBrowser for help on using the repository browser.