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 @ 9332

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

JPALM -- MEDUSA Mocsy output check

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