source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 10047

Last change on this file since 10047 was 10047, checked in by jpalmier, 3 years ago

merge with GO6_package_branch 9385-10020 ; plus debug OMIP_DIC

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