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

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

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

correct PI CO2 flux bug

File size: 39.4 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#   if defined key_omip_dic
460               !! PI DIC - for omip runs
461               f_pi_co2flux(ji,jj) = (1. - fr_i(ji,jj)) * f_pi_co2flux(ji,jj) * &
462                                      86400. / fse3t(ji,jj,1)
463#   endif
464               !! oxygen (O2); OCMIP-2 code
465               !! AXY (23/06/15): amend input list for oxygen to account
466               !!                 for common gas transfer velocity
467               CALL trc_oxy_medusa(ztmp(ji,jj),zsal(ji,jj),f_kw660(ji,jj),     &
468                                   f_pp0(ji,jj),zoxy(ji,jj),                   &
469                                   f_kwo2_dum,f_o2flux(ji,jj),                 &
470                                   f_o2sat(ji,jj))
471               !!
472               !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide
473               !! through by layer thickness
474               f_o2flux(ji,jj)  = (1. - fr_i(ji,jj)) * f_o2flux(ji,jj) *       &
475                                  86400. / fse3t(ji,jj,1)
476            ENDIF
477         ENDDO
478      ENDDO
479
480      !! Jpalm (08-2014)
481      !! DMS surface concentration calculation
482      !! initialy added for UKESM1 model.
483      !! using MET-OFFICE subroutine.
484      !! DMS module only needs Chl concentration and MLD
485      !! to get an aproximate value of DMS concentration.
486      !! air-sea fluxes are calculated by atmospheric chemitry model
487      !! from atm and oc-surface concentrations.
488      !!
489      !! AXY (13/03/15): this is amended to calculate all of the DMS
490      !!                 estimates examined during UKESM1 (see
491      !!                 comments in trcdms_medusa.F90)
492      !!
493      !! AXY (25/05/17): amended to additionally pass DIN limitation as well as [DIN];
494      !!                 accounts for differences in nutrient half-saturations; changes
495      !!                 also made in trc_dms_medusa; this permits an additional DMS
496      !!                 calculation while retaining the existing Anderson one
497      !!
498      IF (jdms == 1) THEN
499         DO jj = 2,jpjm1
500            DO ji = 2,jpim1
501               if (tmask(ji,jj,1) == 1) then
502                  !! calculate weighted half-saturation for DIN uptake
503                  dms_wtkn(ji,jj) = ((zphn(ji,jj) * xnln) +            &
504                                     (zphd(ji,jj) * xnld)) /           &
505                                     (zphn(ji,jj) + zphd(ji,jj))       
506                  !!
507                  !! feed in correct inputs
508                  if (jdms_input .eq. 0) then
509                     !! use instantaneous inputs
510                     dms_nlim(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + dms_wtkn(ji,jj))
511                     !!
512                     CALL trc_dms_medusa(zchn(ji,jj),zchd(ji,jj),             &
513                                         hmld(ji,jj),qsr(ji,jj),              &
514                                         zdin(ji,jj), dms_nlim(ji,jj),        &
515                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
516                                         dms_andm)
517                  else
518                     !! use diel-average inputs
519                     dms_nlim(ji,jj) = zn_dms_din(ji,jj) /                    &
520                                      (zn_dms_din(ji,jj) + dms_wtkn(ji,jj))
521                     !!
522                     CALL trc_dms_medusa(zn_dms_chn(ji,jj),zn_dms_chd(ji,jj), &
523                                         zn_dms_mld(ji,jj),zn_dms_qsr(ji,jj), &
524                                         zn_dms_din(ji,jj),dms_nlim(ji,jj),   &
525                                         dms_andr,dms_simo,dms_aran,dms_hall, & 
526                                         dms_andm)
527                  endif
528                  !!
529                  !! assign correct output to variable passed to atmosphere
530                  if (jdms_model .eq. 1) then
531                     dms_surf = dms_andr
532                  elseif (jdms_model .eq. 2) then
533                     dms_surf = dms_simo
534                  elseif (jdms_model .eq. 3) then
535                     dms_surf = dms_aran
536                  elseif (jdms_model .eq. 4) then
537                     dms_surf = dms_hall
538                  elseif (jdms_model .eq. 5) then
539                     dms_surf = dms_andm
540                  endif
541                  !!
542                  !! 2D diag through iom_use
543                  IF( med_diag%DMS_SURF%dgsave ) THEN
544                     dms_surf2d(ji,jj) = dms_surf
545                  ENDIF
546                  IF( med_diag%DMS_ANDR%dgsave ) THEN
547                     dms_andr2d(ji,jj) = dms_andr
548                  ENDIF
549                  IF( med_diag%DMS_SIMO%dgsave ) THEN
550                     dms_simo2d(ji,jj) = dms_simo
551                  ENDIF
552                  IF( med_diag%DMS_ARAN%dgsave ) THEN
553                     dms_aran2d(ji,jj) = dms_aran
554                  ENDIF
555                  IF( med_diag%DMS_HALL%dgsave ) THEN
556                     dms_hall2d(ji,jj) = dms_hall
557                  ENDIF
558                  IF( med_diag%DMS_ANDM%dgsave ) THEN
559                     dms_andm2d(ji,jj) = dms_andm
560                  ENDIF
561               ENDIF
562            ENDDO
563         ENDDO
564#   if defined key_debug_medusa
565         IF (lwp) write (numout,*)                                &
566            'air-sea: finish calculating dms kt = ',kt
567            CALL flush(numout)
568#   endif
569      ENDIF                  !! End IF (jdms == 1)
570
571      !!
572      !! store 2D outputs
573      !!
574      !! JPALM -- 17-11-16 -- put fgco2 out of diag request
575      !!       is needed for coupling; pass through restart
576      DO jj = 2,jpjm1
577         DO ji = 2,jpim1
578            if (tmask(ji,jj,1) == 1) then
579               !! IF( med_diag%FGCO2%dgsave ) THEN
580               !! convert from  mol/m2/day to kg/m2/s
581               !! mmol-C/m3/d -> kg-CO2/m2/s
582               fgco2(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1) *             &
583                              CO2flux_conv
584               !! ENDIF
585               IF ( lk_iomput ) THEN
586                  IF( med_diag%ATM_PCO2%dgsave ) THEN
587                     f_pco2a2d(ji,jj) = f_pco2atm(ji,jj)
588                  ENDIF
589                  IF( med_diag%OCN_PCO2%dgsave ) THEN
590                     f_pco2w2d(ji,jj) = f_pco2w(ji,jj)
591                  ENDIF
592                  IF( med_diag%CO2FLUX%dgsave ) THEN
593                     !! mmol/m3/d -> mmol/m2/d
594                     f_co2flux2d(ji,jj) = f_co2flux(ji,jj) * fse3t(ji,jj,1)
595                  ENDIF
596                  IF( med_diag%TCO2%dgsave ) THEN
597                     f_TDIC2d(ji,jj) = f_TDIC(ji,jj)
598                  ENDIF
599                  IF( med_diag%TALK%dgsave ) THEN
600                     f_TALK2d(ji,jj) = f_TALK(ji,jj)
601                  ENDIF
602                  IF( med_diag%KW660%dgsave ) THEN
603                      f_kw6602d(ji,jj) = f_kw660(ji,jj)
604                  ENDIF
605                  IF( med_diag%ATM_PP0%dgsave ) THEN
606                      f_pp02d(ji,jj) = f_pp0(ji,jj)
607                  ENDIF
608                  IF( med_diag%O2FLUX%dgsave ) THEN
609                     f_o2flux2d(ji,jj) = f_o2flux(ji,jj)
610                  ENDIF
611                  IF( med_diag%O2SAT%dgsave ) THEN
612                     f_o2sat2d(ji,jj) = f_o2sat(ji,jj)
613                  ENDIF
614                  !! AXY (24/11/16): add in extra MOCSY diagnostics
615                  IF( med_diag%ATM_XCO2%dgsave ) THEN
616                     f_xco2a_2d(ji,jj) = f_xco2a(ji,jj)
617                  ENDIF
618                  IF( med_diag%OCN_FCO2%dgsave ) THEN
619                     f_fco2w_2d(ji,jj) = f_fco2w(ji,jj)
620                  ENDIF
621                  IF( med_diag%ATM_FCO2%dgsave ) THEN
622                     f_fco2a_2d(ji,jj) = f_fco2atm(ji,jj)
623                  ENDIF
624                  IF( med_diag%OCN_RHOSW%dgsave ) THEN
625                     f_ocnrhosw_2d(ji,jj) = f_rhosw(ji,jj)
626                  ENDIF
627                  IF( med_diag%OCN_SCHCO2%dgsave ) THEN
628                     f_ocnschco2_2d(ji,jj) = f_schmidtco2(ji,jj)
629                  ENDIF
630                  IF( med_diag%OCN_KWCO2%dgsave ) THEN
631                     f_ocnkwco2_2d(ji,jj) = f_kwco2(ji,jj)
632                  ENDIF
633                  IF( med_diag%OCN_K0%dgsave ) THEN
634                     f_ocnk0_2d(ji,jj) = f_K0(ji,jj)
635                  ENDIF
636                  IF( med_diag%CO2STARAIR%dgsave ) THEN
637                     f_co2starair_2d(ji,jj) = f_co2starair(ji,jj)
638                  ENDIF
639                  IF( med_diag%OCN_DPCO2%dgsave ) THEN
640                     f_ocndpco2_2d(ji,jj) = f_dpco2(ji,jj)
641                  ENDIF
642#   if defined key_omip_dic
643                  !! AXY (06/08/18): diagnostics for OMIP PI DIC
644                  IF( med_diag%PI_ATM_PCO2%dgsave ) THEN
645                     f_pi_pco2a_2d(ji,jj) = f_pi_pco2atm(ji,jj)
646                  ENDIF
647                  IF( med_diag%PI_OCN_PCO2%dgsave ) THEN
648                     f_pi_pco2w_2d(ji,jj) = f_pi_pco2w(ji,jj)
649                  ENDIF
650                  IF( med_diag%PI_CO2FLUX%dgsave ) THEN
651                     !! mmol/m3/d -> mmol/m2/d
652                     f_pi_co2flux_2d(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1)
653                  ENDIF
654                  IF( med_diag%PI_TCO2%dgsave ) THEN
655                     f_PI_TDIC_2d(ji,jj) = f_PI_TDIC(ji,jj)
656                  ENDIF
657                  IF( med_diag%PI_ATM_XCO2%dgsave ) THEN
658                     f_pi_xco2a_2d(ji,jj) = f_pi_xco2a(ji,jj)
659                  ENDIF
660                  IF( med_diag%PI_OCN_FCO2%dgsave ) THEN
661                     f_pi_fco2w_2d(ji,jj) = f_pi_fco2w(ji,jj)
662                  ENDIF
663                  IF( med_diag%PI_ATM_FCO2%dgsave ) THEN
664                     f_pi_fco2a_2d(ji,jj) = f_pi_fco2atm(ji,jj)
665                  ENDIF
666                  IF( med_diag%PI_CO2STARAIR%dgsave ) THEN
667                     f_pi_co2starair_2d(ji,jj) = f_pi_co2starair(ji,jj)
668                  ENDIF
669                  IF( med_diag%PI_OCN_DPCO2%dgsave ) THEN
670                     f_pi_ocndpco2_2d(ji,jj) = f_pi_dpco2(ji,jj)
671                  ENDIF                 
672                  IF( med_diag%PI_FGCO2%dgsave ) THEN
673                     !! convert from  mol/m2/day to kg/m2/s
674                     !! mmol-C/m3/d -> kg-CO2/m2/s
675                     pi_fgco2(ji,jj) = f_pi_co2flux(ji,jj) * fse3t(ji,jj,1) * &
676                          CO2flux_conv
677                  ENDIF
678                  IF( med_diag%PI_OCN_PH%dgsave ) THEN
679                     f_pi_ph_2d(ji,jj) = f_pi_ph(ji,jj)
680                  ENDIF                 
681                  IF( med_diag%PI_OCNH2CO3%dgsave ) THEN
682                     f_pi_h2co3_2d(ji,jj) = f_pi_h2co3(ji,jj)
683                  ENDIF                 
684                  IF( med_diag%PI_OCN_HCO3%dgsave ) THEN
685                     f_pi_hco3_2d(ji,jj) = f_pi_hco3(ji,jj)
686                  ENDIF                 
687                  IF( med_diag%PI_OCN_CO3%dgsave ) THEN
688                     f_pi_co3_2d(ji,jj) = f_pi_co3(ji,jj)
689                  ENDIF                 
690                  IF( med_diag%PI_OM_CAL%dgsave ) THEN
691                     f_pi_omcal_2d(ji,jj) = f_pi_omcal(ji,jj)
692                  ENDIF                 
693                  IF( med_diag%PI_OM_ARG%dgsave ) THEN
694                     f_pi_omarg_2d(ji,jj) = f_pi_omarg(ji,jj)
695                  ENDIF                 
696#   endif                 
697               ENDIF
698            ENDIF
699         ENDDO
700      ENDDO
701
702      !!-----------------------------------------------------------------
703      !! River inputs
704      !!-----------------------------------------------------------------
705      DO jj = 2,jpjm1
706         DO ji = 2,jpim1
707            !! OPEN wet point IF..THEN loop
708            if (tmask(ji,jj,1) == 1) then
709               !!
710               !! runoff comes in as        kg / m2 / s
711               !! used and written out as   m3 / m2 / d (= m / d)
712               !! where                     1000 kg / m2 / d =
713               !!                             1 m3 / m2 / d = 1 m / d
714               !!
715               !! AXY (17/07/14): the compiler doesn't like this line for
716               !!                 some reason; as MEDUSA doesn't even use
717               !!                 runoff for riverine inputs, a temporary
718               !!                 solution is to switch off runoff entirely
719               !!                 here; again, this change is one of several
720               !!                 that will need revisiting once MEDUSA has
721               !!                 bedded down in UKESM1; particularly so if
722               !!                 the land scheme provides information
723               !!                 concerning nutrient fluxes
724               !!
725               !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. *   &
726               !!                   60. * 24.
727               f_runoff(ji,jj) = 0.0
728               !!
729               !! nutrients are added via rivers to the model in one of
730               !! two ways:
731               !!   1. via river concentration; i.e. the average nutrient
732               !!      concentration of a river water is described by a
733               !!      spatial file, and this is multiplied by runoff to
734               !!      give a nutrient flux
735               !!   2. via direct river flux; i.e. the average nutrient
736               !!      flux due to rivers is described by a spatial file,
737               !!      and this is simply applied as a direct nutrient
738               !!      flux (i.e. it does not relate or respond to model
739               !!      runoff) nutrient fields are derived from the
740               !!      GlobalNEWS 2 database; carbon and alkalinity are
741               !!      derived from continent-scale DIC estimates (Huang et
742               !!      al., 2012) and some Arctic river alkalinity
743               !!      estimates (Katya?)
744               !!
745               !! as of 19/07/12, riverine nutrients can now be spread
746               !! vertically across several grid cells rather than just
747               !! poured into the surface box; this block of code is still
748               !! executed, however, to set up the total amounts of
749               !! nutrient entering via rivers
750               !!
751               !! nitrogen
752               if (jriver_n .eq. 1) then
753                  !! river concentration specified; use runoff to
754                  !! calculate input
755                  f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj)
756               elseif (jriver_n .eq. 2) then
757                  !! river flux specified; independent of runoff
758                  f_riv_n(ji,jj) = riv_n(ji,jj)
759               endif
760               !!
761               !! silicon
762               if (jriver_si .eq. 1) then
763                  !! river concentration specified; use runoff to
764                  !! calculate input
765                  f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj)
766               elseif (jriver_si .eq. 2) then
767                  !! river flux specified; independent of runoff
768                  f_riv_si(ji,jj) = riv_si(ji,jj)
769               endif
770               !!
771               !! carbon
772               if (jriver_c .eq. 1) then
773                  !! river concentration specified; use runoff to
774                  !! calculate input
775                  f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj)
776               elseif (jriver_c .eq. 2) then
777                  !! river flux specified; independent of runoff
778                  f_riv_c(ji,jj) = riv_c(ji,jj)
779               endif
780               !!
781               !! alkalinity
782               if (jriver_alk .eq. 1) then
783                  !! river concentration specified; use runoff to
784                  !! calculate input
785                  f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj)
786               elseif (jriver_alk .eq. 2) then
787                  !! river flux specified; independent of runoff
788                  f_riv_alk(ji,jj) = riv_alk(ji,jj)
789               endif
790            ENDIF
791         ENDDO
792      ENDDO
793# endif  !! Roam
794     
795   END SUBROUTINE air_sea
796
797#else  !! medusa
798   !!======================================================================
799   !!  Dummy module :                                   No MEDUSA bio-model
800   !!======================================================================
801CONTAINS
802   SUBROUTINE air_sea( )                    ! Empty routine
803      WRITE(*,*) 'air_sea: You should not have seen this print! error?'
804   END SUBROUTINE air_sea
805#endif  !! medusa
806
807   !!======================================================================
808END MODULE air_sea_mod
Note: See TracBrowser for help on using the repository browser.