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

source: branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 9073

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

add all micro boil checks and securities

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