source: branches/NERC/dev_r5518_GO6_Carb_Debug/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 8642

Last change on this file since 8642 was 8642, checked in by jpalmier, 4 years ago

jpalm —19-10-17— Carb failure due to 1-cell exceptionnal and ephemeral T increase - update T passed to MEDUSA and add kill switch

File size: 32.0 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 oce,               ONLY: PCO2a_in_cpl
70      USE par_kind,          ONLY: wp
71      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1
72      USE sbc_oce,           ONLY: fr_i, lk_oasis, qsr, wndm
73      USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          &
74                                   jriver_alk, jriver_c,                  &
75                                   jriver_n, jriver_si,                   &
76                                   riv_alk, riv_c, riv_n, riv_si,         &
77                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    &
78                                   zn_dms_mld, zn_dms_qsr,                &
79                                   xnln, xnld 
80      USE trc,               ONLY: med_diag
81      USE zdfmxl,            ONLY: hmld
82
83# if defined key_roam
84      USE gastransfer,       ONLY: gas_transfer
85#  if defined key_mocsy
86      USE mocsy_wrapper,     ONLY: mocsy_interface
87#  else
88      USE trcco2_medusa,     ONLY: trc_co2_medusa
89#  endif
90      USE trcdms_medusa,     ONLY: trc_dms_medusa
91      USE trcoxy_medusa,     ONLY: trc_oxy_medusa
92# endif
93      USE lib_mpp,           ONLY: ctl_stop
94      USE trcstat 
95
96   !!* Substitution
97#  include "domzgr_substitute.h90"
98
99      !! time (integer timestep)
100      INTEGER, INTENT( in ) :: kt
101
102      !! Loop variables
103      INTEGER :: ji, jj
104
105# if defined key_roam
106      !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s
107      REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol
108      REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d
109      REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day
110
111      INTEGER :: iters
112
113      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen
114      REAL(wp), DIMENSION(jpi,jpj) :: f_fco2w, f_rhosw
115      REAL(wp), DIMENSION(jpi,jpj) :: f_fco2atm
116      REAL(wp), DIMENSION(jpi,jpj) :: f_schmidtco2, f_kwco2, f_K0
117      REAL(wp), DIMENSION(jpi,jpj) :: f_co2starair, f_dpco2
118      !! Output arguments from mocsy_interface, which aren't used
119      REAL(wp) :: f_BetaD_dum, f_opres_dum
120      REAL(wp) :: f_insitut_dum
121      REAL(wp) :: f_kwo2_dum
122# endif
123
124
125# if defined key_roam
126      !! init
127      f_fco2w(:,:)       = 0.0
128      f_fco2atm(:,:)     = 0.0
129      f_schmidtco2(:,:)  = 0.0
130      f_kwco2(:,:)       = 0.0
131      f_co2starair(:,:)  = 0.0
132      f_dpco2(:,:)       = 0.0
133      f_rhosw(:,:)       = 0.0
134      f_K0(:,:)          = 0.0
135      !! air pressure (atm); ultimately this will use air
136      !! pressure at the base of the UKESM1 atmosphere
137      !!                                     
138      f_pp0(:,:)   = 1.0
139
140
141      !!-----------------------------------------------------------
142      !! Air-sea gas exchange
143      !!-----------------------------------------------------------
144
145#   if defined key_debug_medusa
146               IF (lwp) write (numout,*)                     & 
147               'air-sea: gas_transfer kt = ', kt
148               CALL flush(numout)
149#   endif
150      DO jj = 2,jpjm1
151         DO ji = 2,jpim1
152            !! OPEN wet point IF..THEN loop
153            IF (tmask(ji,jj,1) == 1) then
154               IF (lk_oasis) THEN
155                  !! use 2D atm xCO2 from atm coupling
156                  f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj)
157                  !!!
158                  !!! Jpalm test on atm xCO2
159                  IF ( (f_xco2a(ji,jj) > 1000 ).OR.(f_xco2a(ji,jj) < 100 ) ) THEN
160                     IF(lwp) WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),   &
161                                             ' -- ji =', mig(ji),' jj = ', mjg(jj)
162                     CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic atm xCO2 ' )
163                 ENDIF
164               ENDIF
165               !!
166               !! AXY (23/06/15): as part of an effort to update the
167               !!                 carbonate chemistry in MEDUSA, the gas
168               !!                 transfer velocity used in the carbon
169               !!                 and oxygen cycles has been harmonised
170               !!                 and is calculated by the same function
171               !!                 here; this harmonisation includes
172               !!                 changes to the PML carbonate chemistry
173               !!                 scheme so that it too makes use of the
174               !!                 same gas transfer velocity; the
175               !!                 preferred parameterisation of this is
176               !!                 Wanninkhof (2014), option 7
177               !!
178               CALL gas_transfer( wndm(ji,jj), 1, 7,         &  ! inputs
179                                  f_kw660(ji,jj) )              ! outputs
180            ENDIF
181         ENDDO
182      ENDDO
183
184#   if defined key_debug_medusa
185               IF (lwp) write (numout,*)                     &
186               'air-sea: carb-chem kt = ', kt
187               CALL flush(numout)
188               !! JPALM add carb print:
189               call trc_rst_dia_stat(f_xco2a(:,:), 'f_xco2a')
190               call trc_rst_dia_stat(wndm(:,:), 'wndm')
191               call trc_rst_dia_stat(f_kw660(:,:), 'f_kw660')
192               call trc_rst_dia_stat(ztmp(:,:), 'ztmp')
193               call trc_rst_dia_stat(zsal(:,:), 'zsal')
194               call trc_rst_dia_stat(zalk(:,:), 'zalk')
195               call trc_rst_dia_stat(zdic(:,:), 'zdic')
196               call trc_rst_dia_stat(zsil(:,:), 'zsil')
197               call trc_rst_dia_stat(zpho(:,:), 'zpho')
198#   endif
199      DO jj = 2,jpjm1
200         DO ji = 2,jpim1
201            if (tmask(ji,jj,1) == 1) then
202               !!
203#  if defined key_axy_carbchem
204#   if defined key_mocsy
205               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
206               !!          carbonate failure in the coupled run.
207               !!          must be associated to air-sea flux or air xCO2...i
208               !!          Check MOCSY inputs
209               IF ( (zsal(ji,jj) > 75.0 ).OR.(zsal(ji,jj) < 0.0 ) .OR.     &
210                    (ztmp(ji,jj) > 50.0 ).OR.(ztmp(ji,jj) < -20.0 ) .OR.     &
211                    (zalk(ji,jj) > 35.0E2 ).OR.(zalk(ji,jj) <= 0.0 ) .OR.     &
212                    (zdic(ji,jj) > 35.0E2 ).OR.(zdic(ji,jj) <= 0.0 ) .OR.     &
213                    (f_kw660(ji,jj) > 1.0E-2 ).OR.(f_kw660(ji,jj) < 0.0 ) ) THEN
214                  IF(lwp) THEN
215                      WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
216                      WRITE(numout,*) ' surface S = ',zsal(ji,jj)
217                      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
218                      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
219                      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
220                      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
221                      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
222                      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
223                      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
224                      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
225                      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
226                      WRITE(numout,*) ' MOCSY input: ji =', mig(ji),' jj = ', mjg(jj),  &
227                                       ' kt = ', kt 
228                      WRITE(numout,*) 'MEDUSA - Air-Sea : unrealistic surface Carb. Chemistry'
229                      CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic surface Carb. Chemistry -- INPUTS' )
230                  ENDIF     
231               ENDIF     
232               !!
233               !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate
234               !!                 chemistry package; note that depth is set to
235               !!                 zero in this call
236               CALL mocsy_interface(ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj),     &
237                                    zdic(ji,jj),zsil(ji,jj),zpho(ji,jj),     &
238                                    f_pp0(ji,jj),0.0,                        &
239                                    gphit(ji,jj),f_kw660(ji,jj),             &
240                                    f_xco2a(ji,jj),1,f_ph(ji,jj),            &
241                                    f_pco2w(ji,jj),f_fco2w(ji,jj),           &
242                                    f_h2co3(ji,jj),f_hco3(ji,jj),            &
243                                    f_co3(ji,jj),f_omarg(ji,jj),             &
244                                    f_omcal(ji,jj),f_BetaD_dum,              &
245                                    f_rhosw(ji,jj),f_opres_dum,              &
246                                    f_insitut_dum,f_pco2atm(ji,jj),          &
247                                    f_fco2atm(ji,jj),f_schmidtco2(ji,jj),    &
248                                    f_kwco2(ji,jj),f_K0(ji,jj),              &
249                                    f_co2starair(ji,jj),f_co2flux(ji,jj),    &
250                                    f_dpco2(ji,jj))
251               !! mmol / m3 -> umol / kg
252               f_TDIC(ji,jj) = (zdic(ji,jj) / f_rhosw(ji,jj)) * 1000.
253               !! meq / m3 ->  ueq / kg
254               f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000.
255               f_dcf(ji,jj)  = f_rhosw(ji,jj)
256               !! Jpalm -- 12-09-2017 -- add extra check after reccurent
257               !!          carbonate failure in the coupled run.
258               !!          must be associated to air-sea flux or air xCO2...i
259               !!          Check MOCSY inputs
260              IF ( (f_pco2w(ji,jj) > 1.E4 ).OR.(f_pco2w(ji,jj) < 0.0 ) .OR.     &
261                   (f_fco2w(ji,jj) > 1.E4 ).OR.(f_fco2w(ji,jj) < 0.0 ) .OR.     &   
262                   (f_fco2atm(ji,jj) > 1.E4 ).OR.(f_fco2atm(ji,jj) < 0.0 ) .OR.     &
263                   (f_co2flux(ji,jj) > 1.E-2 ).OR.(f_co2flux(ji,jj) < -1.E-2 ) .OR.     &
264                   (f_dpco2(ji,jj) > 1.E4 ).OR.(f_dpco2(ji,jj) < -1.E4 ) ) THEN
265                 IF(lwp) THEN
266                     WRITE(numout,*) ' surface T = ',ztmp(ji,jj)
267                     WRITE(numout,*) ' surface S = ',zsal(ji,jj)
268                     WRITE(numout,*) ' surface ALK = ',zalk(ji,jj)
269                     WRITE(numout,*) ' surface DIC = ',zdic(ji,jj)
270                     WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj)
271                     WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)   
272                     WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj)
273                     WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj)
274                     WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj)
275                     WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj)
276                     WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj)
277                     WRITE(numout,*) ' MOCSY output: ji =', mig(ji),' jj = ', mjg(jj),  &
278                                       ' kt = ', kt     
279                     WRITE(numout,*) 'MEDUSA - Air-Sea : unrealistic surface Carb. Chemistry'
280                     CALL ctl_stop( 'MEDUSA - Air-Sea :', 'unrealistic surface Carb. Chemistry -- OUTPUTS' )
281                 ENDIF     
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.