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

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

JPALM — 19-10-17 — Add ctl_warn for master node output

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