source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 10045

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

Andrew's changes to add the OMIP double_DIC (activated with key_omip_dic)

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