source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 10020

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

GMED ticket 406. CPP key fixes.

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