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

Last change on this file since 9257 was 9257, checked in by frrh, 3 years ago

Commit JP's Met Office GMED ticket 371 for trapping or notifying of
peculiar values of MEDUSA fields arising from transient temperature
spikes in the ocean.

Committed using:
svn merge: 9177:9249 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_9163

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