New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
air_sea.F90 in branches/UKMO/dev_r5518_GO6_fix_key_comp/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_fix_key_comp/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 9991

Last change on this file since 9991 was 9991, checked in by frrh, 6 years ago

Fixes to allow MEDUSA to compile with C1D without
the need for multiple (apparently) unrelated CPP keys
merely to satisfy spurious code interdependencies.

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