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

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 15473

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

Use cool skin SST in air-sea CO2 flux calculations.

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