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

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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.

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