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_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 12102

Last change on this file since 12102 was 12102, checked in by dford, 3 years ago

Fix syntax error and relax diagnostic check. See Met Office utils ticket 195.

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