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

source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 13709

Last change on this file since 13709 was 9374, checked in by jpalmier, 6 years ago

JPALM -- 06-03-2018 -- corrections after review suggestions ; see GMED 380 - comment 22

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