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

Last change on this file was 12407, checked in by frrh, 7 months ago

Apply fix for DMS coupling in MEDUSA from Met Office
GMED ticket 514. Merge command:
svn merge -r 12215:12308 svn+ssh://frrh@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/NERC/dev_r5518_GO6_DmsOutFix

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