New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
air_sea.F90 in branches/UKMO/dev_r5518_GO6_package_FOAMv14_Hs_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_Hs_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 15463

Last change on this file since 15463 was 15463, checked in by dford, 8 months ago

CO2 flux based on input significant wave height.

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