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

Last change on this file since 10149 was 10149, checked in by frrh, 23 months ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

svn merge -r 10054:10141 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3

File size: 32.8 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                  IF( med_diag%DMS_SURF%dgsave ) THEN
474                     dms_surf2d(ji,jj) = dms_surf
475                  ENDIF
476                  IF( med_diag%DMS_ANDR%dgsave ) THEN
477                     dms_andr2d(ji,jj) = dms_andr
478                  ENDIF
479                  IF( med_diag%DMS_SIMO%dgsave ) THEN
480                     dms_simo2d(ji,jj) = dms_simo
481                  ENDIF
482                  IF( med_diag%DMS_ARAN%dgsave ) THEN
483                     dms_aran2d(ji,jj) = dms_aran
484                  ENDIF
485                  IF( med_diag%DMS_HALL%dgsave ) THEN
486                     dms_hall2d(ji,jj) = dms_hall
487                  ENDIF
488                  IF( med_diag%DMS_ANDM%dgsave ) THEN
489                     dms_andm2d(ji,jj) = dms_andm
490                  ENDIF
491               ENDIF
492            ENDDO
493         ENDDO
494#   if defined key_debug_medusa
495         IF (lwp) write (numout,*)                                &
496            'air-sea: finish calculating dms kt = ',kt
497            CALL flush(numout)
498#   endif
499      ENDIF                  !! End IF (jdms == 1)
500
501      !!
502      !! store 2D outputs
503      !!
504      !! JPALM -- 17-11-16 -- put fgco2 out of diag request
505      !!       is needed for coupling; pass through restart
506      DO jj = 2,jpjm1
507         DO ji = 2,jpim1
508            if (tmask(ji,jj,1) == 1) then
509               !! IF( med_diag%FGCO2%dgsave ) THEN
510               !! convert from  mol/m2/day to kg/m2/s
511               !! mmol-C/m3/d -> kg-CO2/m2/s
512               fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) *             &
513                              CO2flux_conv
514               !! ENDIF
515               IF ( ln_foam_medusa ) THEN
516                  !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator
517                  f2_pco2w(ji,jj) = f_pco2w(ji,jj)
518                  f2_fco2w(ji,jj) = f_pco2w(ji,jj)
519               ENDIF
520               IF ( lk_iomput ) THEN
521                  IF( med_diag%ATM_PCO2%dgsave ) THEN
522                     f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
523                  ENDIF
524                  IF( med_diag%OCN_PCO2%dgsave ) THEN
525                     f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
526                  ENDIF
527                  IF( med_diag%CO2FLUX%dgsave ) THEN
528                     !! mmol/m3/d -> mmol/m2/d
529                     f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1)
530                  ENDIF
531                  IF( med_diag%TCO2%dgsave ) THEN
532                     f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
533                  ENDIF
534                  IF( med_diag%TALK%dgsave ) THEN
535                     f_TALK2d(ji,jj) = f_TALK(ji,jj)
536                  ENDIF
537                  IF( med_diag%KW660%dgsave ) THEN
538                      f_kw6602d(ji,jj) = f_kw660(ji,jj)
539                  ENDIF
540                  IF( med_diag%ATM_PP0%dgsave ) THEN
541                      f_pp02d(ji,jj) = f_pp0(ji,jj)
542                  ENDIF
543                  IF( med_diag%O2FLUX%dgsave ) THEN
544                     f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
545                  ENDIF
546                  IF( med_diag%O2SAT%dgsave ) THEN
547                     f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
548                  ENDIF
549                  !! AXY (24/11/16): add in extra MOCSY diagnostics
550                  IF( med_diag%ATM_XCO2%dgsave ) THEN
551                     f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
552                  ENDIF
553                  IF( med_diag%OCN_FCO2%dgsave ) THEN
554                     f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
555                  ENDIF
556                  IF( med_diag%ATM_FCO2%dgsave ) THEN
557                     f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
558                  ENDIF
559                  IF( med_diag%OCN_RHOSW%dgsave ) THEN
560                     f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
561                  ENDIF
562                  IF( med_diag%OCN_SCHCO2%dgsave ) THEN
563                     f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
564                  ENDIF
565                  IF( med_diag%OCN_KWCO2%dgsave ) THEN
566                     f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
567                  ENDIF
568                  IF( med_diag%OCN_K0%dgsave ) THEN
569                     f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
570                  ENDIF
571                  IF( med_diag%CO2STARAIR%dgsave ) THEN
572                     f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
573                  ENDIF
574                  IF( med_diag%OCN_DPCO2%dgsave ) THEN
575                     f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
576                  ENDIF
577               ENDIF
578            ENDIF
579         ENDDO
580      ENDDO
581
582# endif
583
584      !!-----------------------------------------------------------------
585      !! River inputs
586      !!-----------------------------------------------------------------
587      DO jj = 2,jpjm1
588         DO ji = 2,jpim1
589            !! OPEN wet point IF..THEN loop
590            if (tmask(ji,jj,1) == 1) then
591               !!
592               !! runoff comes in as        kg / m2 / s
593               !! used and written out as   m3 / m2 / d (= m / d)
594               !! where                     1000 kg / m2 / d =
595               !!                             1 m3 / m2 / d = 1 m / d
596               !!
597               !! AXY (17/07/14): the compiler doesn't like this line for
598               !!                 some reason; as MEDUSA doesn't even use
599               !!                 runoff for riverine inputs, a temporary
600               !!                 solution is to switch off runoff entirely
601               !!                 here; again, this change is one of several
602               !!                 that will need revisiting once MEDUSA has
603               !!                 bedded down in UKESM1; particularly so if
604               !!                 the land scheme provides information
605               !!                 concerning nutrient fluxes
606               !!
607               !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. *   &
608               !!                   60. * 24.
609               f_runoff(ji,jj) = 0.0
610               !!
611               !! nutrients are added via rivers to the model in one of
612               !! two ways:
613               !!   1. via river concentration; i.e. the average nutrient
614               !!      concentration of a river water is described by a
615               !!      spatial file, and this is multiplied by runoff to
616               !!      give a nutrient flux
617               !!   2. via direct river flux; i.e. the average nutrient
618               !!      flux due to rivers is described by a spatial file,
619               !!      and this is simply applied as a direct nutrient
620               !!      flux (i.e. it does not relate or respond to model
621               !!      runoff) nutrient fields are derived from the
622               !!      GlobalNEWS 2 database; carbon and alkalinity are
623               !!      derived from continent-scale DIC estimates (Huang et
624               !!      al., 2012) and some Arctic river alkalinity
625               !!      estimates (Katya?)
626               !!
627               !! as of 19/07/12, riverine nutrients can now be spread
628               !! vertically across several grid cells rather than just
629               !! poured into the surface box; this block of code is still
630               !! executed, however, to set up the total amounts of
631               !! nutrient entering via rivers
632               !!
633               !! nitrogen
634               if (jriver_n .eq. 1) then
635                  !! river concentration specified; use runoff to
636                  !! calculate input
637                  f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
638               elseif (jriver_n .eq. 2) then
639                  !! river flux specified; independent of runoff
640                  f_riv_n(ji,jj) = riv_n(ji,jj)
641               endif
642               !!
643               !! silicon
644               if (jriver_si .eq. 1) then
645                  !! river concentration specified; use runoff to
646                  !! calculate input
647                  f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
648               elseif (jriver_si .eq. 2) then
649                  !! river flux specified; independent of runoff
650                  f_riv_si(ji,jj) = riv_si(ji,jj)
651               endif
652               !!
653               !! carbon
654               if (jriver_c .eq. 1) then
655                  !! river concentration specified; use runoff to
656                  !! calculate input
657                  f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
658               elseif (jriver_c .eq. 2) then
659                  !! river flux specified; independent of runoff
660                  f_riv_c(ji,jj) = riv_c(ji,jj)
661               endif
662               !!
663               !! alkalinity
664               if (jriver_alk .eq. 1) then
665                  !! river concentration specified; use runoff to
666                  !! calculate input
667                  f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
668               elseif (jriver_alk .eq. 2) then
669                  !! river flux specified; independent of runoff
670                  f_riv_alk(ji,jj) = riv_alk(ji,jj)
671               endif
672            ENDIF
673         ENDDO
674      ENDDO
675
676   END SUBROUTINE air_sea
677
678#else
679   !!======================================================================
680   !!  Dummy module :                                   No MEDUSA bio-model
681   !!======================================================================
682CONTAINS
683   SUBROUTINE air_sea( )                    ! Empty routine
684      WRITE(*,*) 'air_sea: You should not have seen this print! error?'
685   END SUBROUTINE air_sea
686#endif 
687
688   !!======================================================================
689END MODULE air_sea_mod
Note: See TracBrowser for help on using the repository browser.