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

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 9292

Last change on this file since 9292 was 9292, checked in by dford, 6 years ago

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc_v2.

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