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.
Changeset 5841 – NEMO

Changeset 5841


Ignore:
Timestamp:
2015-10-30T12:48:06+01:00 (8 years ago)
Author:
jpalmier
Message:

JPALM --30-10-2015-- Add MOCSY and DMS to MEDUSA-NEMO3.6

Location:
branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC
Files:
8 added
14 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r5726 r5841  
    172172!! 
    173173!! UKESM diagnostics 
    174    INTEGER  ::  jdms      !:  include DMS diagnostics ? Jpalm (27-08-2014)  
     174   INTEGER  ::  jdms         !: include DMS diagnostics ? Jpalm (27-08-2014)  
     175   INTEGER  ::  jdms_input   !: use instant (0) or diel-average (1) inputs (AXY, 08/07/2015) 
     176   INTEGER  ::  jdms_model   !: choice of DMS model passed to atmosphere 
     177!!                              1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL 
    175178!! 
    176179!! 
     
    217220   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_sed_ca   !: 2D inorganic carbon   (now) 
    218221   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_sed_ca   !: 2D inorganic carbon   (after) 
     222!! 
     223!! 2D fields of temporally averaged properties for DMS calculations (AXY, 07/07/15) 
     224   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_dms_chn  !: 2D avg CHN   (before) 
     225   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_dms_chn  !: 2D avg CHN   (now) 
     226   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_dms_chn  !: 2D avg CHN   (after) 
     227   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_dms_chd  !: 2D avg CHD   (before) 
     228   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_dms_chd  !: 2D avg CHD   (now) 
     229   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_dms_chd  !: 2D avg CHD   (after) 
     230   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_dms_mld  !: 2D avg MLD   (before) 
     231   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_dms_mld  !: 2D avg MLD   (now) 
     232   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_dms_mld  !: 2D avg MLD   (after) 
     233   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_dms_qsr  !: 2D avg QSR   (before) 
     234   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_dms_qsr  !: 2D avg QSR   (now) 
     235   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_dms_qsr  !: 2D avg QSR   (after) 
     236   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_dms_din  !: 2D avg DIN   (before) 
     237   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zn_dms_din  !: 2D avg DIN   (now) 
     238   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: za_dms_din  !: 2D avg DIN   (after) 
    219239#endif 
    220240 
     
    415435      !!---------------------------------------------------------------------- 
    416436      USE lib_mpp , ONLY: ctl_warn 
    417       INTEGER ::   ierr(6)        ! Local variables 
     437      INTEGER ::   ierr(7)        ! Local variables 
    418438      !!---------------------------------------------------------------------- 
    419439      ierr(:) = 0 
     
    439459         &      zb_sed_ca(jpi,jpj)   , zn_sed_ca(jpi,jpj)   ,       & 
    440460         &      za_sed_ca(jpi,jpj)   ,                           STAT=ierr(3) ) 
     461      !* 2D fields of temporally averaged properties for DMS calculations (AXY, 07/07/15) 
     462      ALLOCATE( zb_dms_chn(jpi,jpj)  , zn_dms_chn(jpi,jpj)  ,       & 
     463         &      za_dms_chn(jpi,jpj)  ,                              & 
     464         &      zb_dms_chd(jpi,jpj)  , zn_dms_chd(jpi,jpj)  ,       &         
     465         &      za_dms_chd(jpi,jpj)  ,                              & 
     466         &      zb_dms_mld(jpi,jpj)  , zn_dms_mld(jpi,jpj)  ,       &         
     467         &      za_dms_mld(jpi,jpj)  ,                              & 
     468         &      zb_dms_qsr(jpi,jpj)  , zn_dms_qsr(jpi,jpj)  ,       &         
     469         &      za_dms_qsr(jpi,jpj)  ,                              & 
     470         &      zb_dms_din(jpi,jpj)  , zn_dms_din(jpi,jpj)  ,       &         
     471         &      za_dms_din(jpi,jpj)  ,                           STAT=ierr(4) ) 
    441472# endif 
    442473      !* 2D fields of miscellaneous parameters 
     
    444475         &      dustmo(jpi,jpj,2)    , riv_n(jpi,jpj)       ,       & 
    445476         &      riv_si(jpi,jpj)      , riv_c(jpi,jpj)       ,       & 
    446          &      riv_alk(jpi,jpj)     , friver_dep(jpk,jpk)  ,    STAT=ierr(4) ) 
     477         &      riv_alk(jpi,jpj)     , friver_dep(jpk,jpk)  ,    STAT=ierr(5) ) 
    447478      !* 2D and 3D fields of light parameters 
    448479      ALLOCATE( neln(jpi,jpj)        , xze(jpi,jpj)         ,       & 
    449          &      xpar(jpi,jpj,jpk)    ,                           STAT=ierr(5) ) 
     480         &      xpar(jpi,jpj,jpk)    ,                           STAT=ierr(6) ) 
    450481      !* 2D and 3D fields of sediment-associated parameters 
    451482      ALLOCATE( dminl(jpi,jpj)       , dmin3(jpi,jpj,jpk)   ,       & 
     
    454485         &      fbodf(jpi,jpj)       , fbods(jpi,jpj)       ,       & 
    455486         &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       & 
    456          &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(6) )  
     487         &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(7) )  
    457488#endif 
    458489      ! 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r5726 r5841  
    1616   !!  -   !  2013-05  (A. Yool)              updated for v3.5 
    1717   !!  -   !  2014-08  (A. Yool, J. Palm)     Add DMS module for UKESM1 model 
     18   !!  -   !  2015-06  (A. Yool)              Update to include MOCSY 
     19   !!  -   !  2015-07  (A. Yool)              Update for rolling averages 
    1820   !!---------------------------------------------------------------------- 
    1921   !! 
     
    3638#endif 
    3739   !! 
     40#if defined key_mocsy 
     41   !!---------------------------------------------------------------------- 
     42   !! Updates with the addition of MOCSY include: 
     43   !!   - option to use PML or MOCSY carbonate chemistry (the latter is  
     44   !!     preferred) 
     45   !!   - central calculation of gas transfer velocity, f_kw660; previously 
     46   !!     this was done separately for CO2 and O2 with predictable results 
     47   !!   - distribution of f_kw660 to both PML and MOCSY CO2 air-sea flux  
     48   !!     calculations and to those for O2 air-sea flux 
     49   !!   - extra diagnostics included for MOCSY 
     50   !!---------------------------------------------------------------------- 
     51#endif 
     52   !! 
    3853#if defined key_medusa 
    3954   !!---------------------------------------------------------------------- 
     
    5570# endif 
    5671# if defined key_roam 
     72      USE gastransfer 
     73#  if defined key_mocsy 
     74      USE mocsy_wrapper 
     75#  else 
    5776      USE trcco2_medusa 
     77#  endif 
    5878      USE trcoxy_medusa 
    5979      !! Jpalm (08/08/2014) 
     
    134154      REAL(wp) ::    ztmp, zsal 
    135155# endif 
     156# if defined key_mocsy 
     157      REAL(wp) ::    zpho 
     158# endif 
    136159      !! 
    137160      !! integrated source and sink terms 
     
    266289      REAL(wp) ::    f_kw660, f_o2flux, f_o2sat 
    267290      REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg 
     291      !! 
     292      !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen 
     293      REAL(wp) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm 
     294      REAL(wp) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 
     295      !! 
    268296      INTEGER  ::    iters 
    269297      REAL(wp) ::    f_year 
     
    305333 
    306334      !!--------------------------------------------------------------------- 
     335 
     336# if defined key_debug_medusa 
     337      IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined' 
     338      CALL flush(numout) 
     339# endif  
    307340 
    308341      !! AXY (20/11/14): alter this to report on first MEDUSA call 
     
    486519# endif 
    487520 
     521# if defined key_debug_medusa 
     522      IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 
     523      CALL flush(numout) 
     524# endif  
     525 
    488526# if defined key_roam 
    489527      !!---------------------------------------------------------------------- 
     
    514552         f_pco2a = fq4 
    515553      endif 
    516 # if defined key_axy_pi_co2 
     554#  if defined key_axy_pi_co2 
    517555      f_pco2a = hist_pco2(1) 
    518       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  = FIXED' 
    519 # endif 
     556#  endif 
    520557      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    521558      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day) 
     
    529566# endif 
    530567 
     568# if defined key_debug_medusa 
     569      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry' 
     570      IF (lwp) write (numout,*) 'trc_bio_medusa: kt = ', kt 
     571      IF (lwp) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000 
     572      CALL flush(numout) 
     573# endif  
     574 
    531575# if defined key_roam 
    532       !! AXY (20/11/14): alter to call on first MEDUSA timestep 
     576      !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every 
     577      !!                 month (this is hardwired as 960 timesteps but should 
     578      !!                 be calculated and done properly 
    533579      !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN 
    534       IF( kt == nittrc000 .or. mod(kt,1920) == 0 ) THEN 
     580      IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN 
    535581         !!---------------------------------------------------------------------- 
    536582         !! Calculate the carbonate chemistry for the whole ocean on the first 
     
    540586         !! 
    541587         IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt 
     588         CALL flush(numout) 
    542589         !! blank flags 
    543590         i2_omcal(:,:) = 0 
     
    549596                  !! OPEN wet point IF..THEN loop 
    550597                  if (tmask(ji,jj,jk).eq.1) then 
    551                      !! carbonate chemistry 
     598                     !! do carbonate chemistry 
     599                     !! 
    552600                     fdep2 = fsdept(ji,jj,jk)           !! set up level midpoint 
     601                     !! 
    553602                     !! set up required state variables 
    554603                     zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
     
    556605                     ztmp = tsn(ji,jj,jk,jp_tem)        !! temperature 
    557606                     zsal = tsn(ji,jj,jk,jp_sal)        !! salinity 
     607#  if defined key_mocsy 
     608                     zsil = max(0.,trn(ji,jj,jk,jpsil))        !! silicic acid 
     609                     zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
     610#  endif 
    558611           !! 
    559612           !! AXY (28/02/14): check input fields 
     
    571624                        ji, ',', jj, ',', jk, ') at time', kt 
    572625                     endif 
     626                     !! 
     627                     !! blank input variables not used at this stage (they relate to air-sea flux) 
     628                     f_kw660 = 1.0 
     629                     f_pp0   = 1.0 
     630                     !! 
    573631                     !! calculate carbonate chemistry at grid cell midpoint 
    574                      CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, 5.0, f_pco2a, &    ! inputs 
    575                      f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),            &    ! outputs 
     632#  if defined key_mocsy 
     633                     !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
     634                     !!                 chemistry package 
     635                     CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,         &    ! inputs 
     636                     f_pp0, fdep2, flatx, f_kw660, f_pco2a, 1,                         &    ! inputs 
     637                     f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),   &    ! outputs 
     638                     f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,             &    ! outputs 
     639                     f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,                &    ! outputs 
     640                     f_co2starair, f_co2flux, f_dpco2 )                                     ! outputs 
     641                     !! 
     642                     f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 
     643                     f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg 
     644                     f_dcf  = f_rhosw 
     645#  else 
     646                     !! AXY (22/06/15): use old PML carbonate chemistry package (the 
     647                     !!                 MEDUSA-2 default) 
     648                     CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660,      &    ! inputs 
     649                     f_pco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),   &    ! outputs 
    576650                     f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters)      ! outputs 
    577651                     !!  
     
    581655                        iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
    582656                     endif 
     657#  endif 
     658                     !! 
    583659                     !! store 3D outputs 
    584660                     f3_pH(ji,jj,jk)    = f_ph 
     
    588664                     f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 
    589665                     f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 
     666                     !! 
    590667                     !! CCD calculation: calcite 
    591668                     if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then 
     
    607684                        i2_omcal(ji,jj)   = 1 
    608685                     endif 
     686                     !! 
    609687                     !! CCD calculation: aragonite 
    610688                     if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then 
     
    633711# endif 
    634712 
     713# if defined key_debug_medusa 
     714      IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations' 
     715      CALL flush(numout) 
     716# endif  
     717 
    635718      !!---------------------------------------------------------------------- 
    636719      !! MEDUSA has unified equation through the water column 
     
    656739            !! OPEN wet point IF..THEN loop 
    657740            if (tmask(ji,jj,jk).eq.1) then                
    658  
    659741               !!====================================================================== 
    660742               !! SETUP LOCAL GRID CELL 
     
    814896                  !!---------------------------------------------------------------------- 
    815897                  !! 
    816                   !! a bit of set up ... 
    817                   ! f_uwind = zwnd_i(ji,jj) 
    818                   ! f_vwind = zwnd_j(ji,jj) 
    819898                  !! AXY (17/07/14): zwind_i and zwind_j do not exist in this 
    820899                  !!                 version of NEMO because it does not include 
     
    827906                  !!                 revisited when MEDUSA properly interacts 
    828907                  !!                 with UKESM1 physics 
    829                   ! f_uwind = zwind_i(ji,jj) 
    830                   ! f_vwind = zwind_j(ji,jj) 
    831                   ! f_wind  = ((f_uwind**2.0) + (f_vwind**2.0))**0.5 
     908                  !! 
    832909                  f_wind  = wndm(ji,jj) 
    833                   !! AXY (17/07/14): the current oxygen code takes in separate 
    834                   !!                 U and V components of the wind; to avoid 
    835                   !!                 the need to change this, calculate these 
    836                   !!                 components based on wndm; again, this is 
    837                   !!                 not ideal, but should suffice as a 
    838                   !!                 temporary measure; in the long-term all 
    839                   !!                 of MEDUSA's air-sea gas exchange terms 
    840                   !!                 will be revisited to ensure that they are 
    841                   !!                 valid and self-consistent; and the CO2 
    842                   !!                 code will be wholly replaced with a more 
    843                   !!                 up-to-date parameterisation 
    844                   f_uwind = ((f_wind**2.0) / 2.0)**0.5 
    845                   f_vwind = f_uwind 
     910                  !! 
     911                  !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
     912                  !!                 in MEDUSA, the gas transfer velocity used in the carbon 
     913                  !!                 and oxygen cycles has been harmonised and is calculated 
     914                  !!                 by the same function here; this harmonisation includes 
     915                  !!                 changes to the PML carbonate chemistry scheme so that 
     916                  !!                 it too makes use of the same gas transfer velocity; the 
     917                  !!                 preferred parameterisation of this is Wanninkhof (2014), 
     918                  !!                 option 7 
     919                  !! 
     920#   if defined key_debug_medusa 
     921                     IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
     922                     CALL flush(numout) 
     923#   endif 
     924                  CALL gas_transfer( f_wind, 1, 7, &  ! inputs 
     925                                     f_kw660 )        ! outputs 
     926#   if defined key_debug_medusa 
     927                     IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
     928                     CALL flush(numout) 
     929#   endif 
     930                  !! 
     931                  !! air pressure (atm); ultimately this will use air pressure at the base 
     932                  !! of the UKESM1 atmosphere  
     933                  !!                                      
    846934                  f_pp0   = 1.0 
     935                  !! 
    847936                  !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp 
    848937                  !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
     
    852941                  !! 
    853942#  if defined key_axy_carbchem 
     943#   if defined key_mocsy 
     944                  !! 
     945                  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
     946                  !!                 chemistry package; note that depth is set to 
     947                  !!                 zero in this call 
     948                  CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,        &  ! inputs 
     949                  f_pp0, 0.0, flatx, f_kw660, f_pco2a, 1,                          &  ! inputs 
     950                  f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),  &  ! outputs 
     951                  f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,            &  ! outputs 
     952                  f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,               &  ! outputs 
     953                  f_co2starair, f_co2flux, f_dpco2 )                                  ! outputs 
     954                  !! 
     955                  f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 
     956                  f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg 
     957                  f_dcf  = f_rhosw 
     958#   else                   
    854959                  iters = 0 
    855960                  !! 
    856961                  !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 
    857                   CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_wind, f_pco2a, &     ! inputs 
    858                   f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),             &     ! outputs 
    859                   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters )       ! outputs 
     962                  CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_pco2a,  &  ! inputs 
     963                  f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),               &  ! outputs 
     964                  f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters )      ! outputs 
    860965                  !! 
    861966                  !! AXY (09/01/14): removed iteration and NaN checks; these have 
     
    869974                     iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
    870975                  endif 
     976#   endif 
    871977#  else 
    872978                  !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
     
    882988                  f_TDIC         = zdic 
    883989                  f_TALK         = zalk 
    884                   f_dcf          = 1. 
     990                  f_dcf          = 1.026 
    885991                  f_henry        = 1. 
     992                  !! AXY (23/06/15): add in some extra MOCSY diagnostics 
     993                  f_fco2w        = fpco2a 
     994                  f_BetaD        = 1. 
     995                  f_rhosw        = 1.026 
     996                  f_opres        = 0. 
     997                  f_insitut      = ztmp 
     998                  f_pco2atm      = fpco2a 
     999                  f_fco2atm      = fpco2a 
     1000                  f_schmidtco2   = 660. 
     1001                  f_kwco2        = 0. 
     1002                  f_K0           = 0. 
     1003                  f_co2starair   = fpco2a 
     1004                  f_dpco2        = 0. 
    8861005#  endif 
    8871006                  !! 
    888                   !! already in right units; correct for sea-ice; divide through by layer thickness 
    889                   f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux / fthk 
     1007                  !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness 
     1008                  f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk 
    8901009                  !! 
    8911010                  !! oxygen (O2); OCMIP-2 code 
    892                   CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, &  ! inputs 
    893                   f_kw660, f_o2flux, f_o2sat )                                                     ! outputs 
    894                   !! 
    895                   !! mol/m3/s -> mmol/m3/d; correct for sea-ice 
    896                   f_o2flux  = (1. - fr_i(ji,jj)) * f_o2flux * 1000. * 60. * 60. * 24.  
    897                   f_o2sat   = f_o2sat  * 1000. 
     1011                  !! AXY (23/06/15): amend input list for oxygen to account for common gas 
     1012                  !!                 transfer velocity 
     1013                  !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk,  &  ! inputs 
     1014                  !! f_kw660, f_o2flux, f_o2sat )                                                      ! outputs 
     1015                  CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy,  &  ! inputs 
     1016                  f_kwo2, f_o2flux, f_o2sat )                                ! outputs 
     1017                  !! 
     1018                  !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness 
     1019                  f_o2flux  = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk 
    8981020                  !! 
    8991021                  !! Jpalm (08-2014) 
     
    9111033                  !! 
    9121034                  IF (jdms .eq. 1) THEN 
    913                   !! 
    914                      !! CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), &                            !! inputs 
    915                      !! dms_surf )                                                                 !! outputs 
    916                      CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, &             !! inputs 
    917                      dms_surf, dms_andr, dms_simo, dms_aran, dms_hall )                            !! outputs 
     1035                     !! 
     1036                     !! feed in correct inputs 
     1037                     if (jdms_input .eq. 0) then 
     1038                        !! use instantaneous inputs 
     1039                        CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, &  ! inputs 
     1040                        dms_andr, dms_simo, dms_aran, dms_hall )                           ! outputs 
     1041                     else 
     1042                        !! use diel-average inputs 
     1043                        CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), &  ! inputs 
     1044                        zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj),   &  ! inputs 
     1045                        dms_andr, dms_simo, dms_aran, dms_hall )                      ! outputs 
     1046                     endif 
     1047                     !! 
     1048                     !! assign correct output to variable passed to atmosphere 
     1049                     if     (jdms_model .eq. 1) then 
     1050                        dms_surf = dms_andr 
     1051                     elseif (jdms_model .eq. 2) then 
     1052                        dms_surf = dms_simo 
     1053                     elseif (jdms_model .eq. 3) then 
     1054                        dms_surf = dms_aran 
     1055                     elseif (jdms_model .eq. 4) then 
     1056                        dms_surf = dms_hall 
     1057                     endif 
    9181058                  ENDIF 
    9191059                  !! End DMS Loop 
     
    9531093                  !! alkalinity are derived from continent-scale DIC estimates (Huang et al.,  
    9541094                  !! 2012) and some Arctic river alkalinity estimates (Katya?) 
    955                   !!       
     1095                  !!  
    9561096                  !! as of 19/07/12, riverine nutrients can now be spread vertically across  
    9571097                  !! several grid cells rather than just poured into the surface box; this 
     
    33173457         ENDDO 
    33183458      ENDDO 
     3459      !! AXY (07/07/15): temporary hijacking 
     3460# if defined key_roam 
     3461      trc2d(:,:,126) = zn_dms_chn(:,:) 
     3462      trc2d(:,:,127) = zn_dms_chd(:,:) 
     3463      trc2d(:,:,128) = zn_dms_mld(:,:) 
     3464      trc2d(:,:,129) = zn_dms_qsr(:,:) 
     3465      trc2d(:,:,130) = zn_dms_din(:,:) 
     3466# endif 
    33193467 
    33203468      if (ibenthic.eq.2) then 
     
    33223470         !! effectively commented out because it does not work as  
    33233471         !! anticipated; it can be deleted at a later date 
    3324       if (jorgben.eq.1) then 
    3325          za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ) * rdt 
    3326          za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt 
    3327          za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ) * rdt 
    3328       endif 
    3329       if (jinorgben.eq.1) then 
    3330          za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt 
    3331          za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt 
    3332       endif 
    3333       !! 
    3334       !! Leap-frog scheme - only in explicit case, otherwise the time stepping 
    3335       !! is already being done in trczdf 
    3336       !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN 
    3337       !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) 
    3338       !!    IF( neuler == 0 .AND. kt == nittrc000 )   zfact = rdttra(jk) * FLOAT(ndttrc) 
    3339       !!    if (jorgben.eq.1) then 
    3340       !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  ) 
    3341       !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) 
    3342       !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  ) 
    3343       !!    endif 
    3344       !!    if (jinorgben.eq.1) then 
    3345       !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) 
    3346       !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) 
    3347       !!    endif 
    3348       !! ENDIF 
    3349       !!  
    3350       !! Time filter and swap of arrays 
    3351       IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN         ! centred or tvd scheme 
    3352          IF( neuler == 0 .AND. kt == nittrc000 ) THEN 
     3472         if (jorgben.eq.1) then 
     3473            za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  ) * rdt 
     3474            za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) ) * rdt 
     3475            za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  ) * rdt 
     3476         endif 
     3477         if (jinorgben.eq.1) then 
     3478            za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt 
     3479            za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt 
     3480         endif 
     3481         !! 
     3482         !! Leap-frog scheme - only in explicit case, otherwise the time stepping 
     3483         !! is already being done in trczdf 
     3484         !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN 
     3485         !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) 
     3486         !!    IF( neuler == 0 .AND. kt == nittrc000 )   zfact = rdttra(jk) * FLOAT(ndttrc) 
     3487         !!    if (jorgben.eq.1) then 
     3488         !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  ) 
     3489         !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) 
     3490         !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  ) 
     3491         !!    endif 
     3492         !!    if (jinorgben.eq.1) then 
     3493         !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) 
     3494         !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) 
     3495         !!    endif 
     3496         !! ENDIF 
     3497         !!  
     3498         !! Time filter and swap of arrays 
     3499         IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme 
     3500            IF( neuler == 0 .AND. kt == nittrc000 ) THEN 
     3501               if (jorgben.eq.1) then 
     3502                  zb_sed_n(:,:)  = zn_sed_n(:,:) 
     3503                  zn_sed_n(:,:)  = za_sed_n(:,:) 
     3504                  za_sed_n(:,:)  = 0.0 
     3505                  !! 
     3506                  zb_sed_fe(:,:) = zn_sed_fe(:,:) 
     3507                  zn_sed_fe(:,:) = za_sed_fe(:,:) 
     3508                  za_sed_fe(:,:) = 0.0 
     3509                  !! 
     3510                  zb_sed_c(:,:)  = zn_sed_c(:,:) 
     3511                  zn_sed_c(:,:)  = za_sed_c(:,:) 
     3512                  za_sed_c(:,:)  = 0.0 
     3513               endif 
     3514               if (jinorgben.eq.1) then 
     3515                  zb_sed_si(:,:) = zn_sed_si(:,:) 
     3516                  zn_sed_si(:,:) = za_sed_si(:,:) 
     3517                  za_sed_si(:,:) = 0.0 
     3518                  !! 
     3519                  zb_sed_ca(:,:) = zn_sed_ca(:,:) 
     3520                  zn_sed_ca(:,:) = za_sed_ca(:,:) 
     3521                  za_sed_ca(:,:) = 0.0 
     3522               endif 
     3523            ELSE 
     3524               if (jorgben.eq.1) then 
     3525                  zb_sed_n(:,:)  = (atfp  * ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) + (atfp1 * zn_sed_n(:,:) ) 
     3526                  zn_sed_n(:,:)  = za_sed_n(:,:) 
     3527                  za_sed_n(:,:)  = 0.0 
     3528                  !! 
     3529                  zb_sed_fe(:,:) = (atfp  * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:)) 
     3530                  zn_sed_fe(:,:) = za_sed_fe(:,:) 
     3531                  za_sed_fe(:,:) = 0.0 
     3532                  !! 
     3533                  zb_sed_c(:,:)  = (atfp  * ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) + (atfp1 * zn_sed_c(:,:) ) 
     3534                  zn_sed_c(:,:)  = za_sed_c(:,:) 
     3535                  za_sed_c(:,:)  = 0.0 
     3536               endif 
     3537               if (jinorgben.eq.1) then 
     3538                  zb_sed_si(:,:) = (atfp  * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:)) 
     3539                  zn_sed_si(:,:) = za_sed_si(:,:) 
     3540                  za_sed_si(:,:) = 0.0 
     3541                  !! 
     3542                  zb_sed_ca(:,:) = (atfp  * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:)) 
     3543                  zn_sed_ca(:,:) = za_sed_ca(:,:) 
     3544                  za_sed_ca(:,:) = 0.0 
     3545               endif 
     3546            ENDIF 
     3547         ELSE                   !  case of smolar scheme or muscl 
    33533548            if (jorgben.eq.1) then 
    3354                zb_sed_n(:,:)  = zn_sed_n(:,:) 
     3549               zb_sed_n(:,:)  = za_sed_n(:,:) 
    33553550               zn_sed_n(:,:)  = za_sed_n(:,:) 
    33563551               za_sed_n(:,:)  = 0.0 
    33573552               !! 
    3358                zb_sed_fe(:,:) = zn_sed_fe(:,:) 
     3553               zb_sed_fe(:,:) = za_sed_fe(:,:) 
    33593554               zn_sed_fe(:,:) = za_sed_fe(:,:) 
    33603555               za_sed_fe(:,:) = 0.0 
    33613556               !! 
    3362                zb_sed_c(:,:)  = zn_sed_c(:,:) 
     3557               zb_sed_c(:,:)  = za_sed_c(:,:) 
    33633558               zn_sed_c(:,:)  = za_sed_c(:,:) 
    33643559               za_sed_c(:,:)  = 0.0 
    33653560            endif 
    33663561            if (jinorgben.eq.1) then 
    3367                zb_sed_si(:,:) = zn_sed_si(:,:) 
     3562               zb_sed_si(:,:) = za_sed_si(:,:) 
    33683563               zn_sed_si(:,:) = za_sed_si(:,:) 
    33693564               za_sed_si(:,:) = 0.0 
    33703565               !! 
    3371                zb_sed_ca(:,:) = zn_sed_ca(:,:) 
    3372                zn_sed_ca(:,:) = za_sed_ca(:,:) 
    3373                za_sed_ca(:,:) = 0.0 
    3374             endif 
    3375          ELSE 
    3376             if (jorgben.eq.1) then 
    3377                zb_sed_n(:,:)  = (atfp  * ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) + (atfp1 * zn_sed_n(:,:) ) 
    3378                zn_sed_n(:,:)  = za_sed_n(:,:) 
    3379                za_sed_n(:,:)  = 0.0 
    3380                !! 
    3381                zb_sed_fe(:,:) = (atfp  * ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) + (atfp1 * zn_sed_fe(:,:)) 
    3382                zn_sed_fe(:,:) = za_sed_fe(:,:) 
    3383                za_sed_fe(:,:) = 0.0 
    3384                !! 
    3385                zb_sed_c(:,:)  = (atfp  * ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) + (atfp1 * zn_sed_c(:,:) ) 
    3386                zn_sed_c(:,:)  = za_sed_c(:,:) 
    3387                za_sed_c(:,:)  = 0.0 
    3388             endif 
    3389             if (jinorgben.eq.1) then 
    3390                zb_sed_si(:,:) = (atfp  * ( zb_sed_si(:,:) + za_sed_si(:,:) )) + (atfp1 * zn_sed_si(:,:)) 
    3391                zn_sed_si(:,:) = za_sed_si(:,:) 
    3392                za_sed_si(:,:) = 0.0 
    3393                !! 
    3394                zb_sed_ca(:,:) = (atfp  * ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) + (atfp1 * zn_sed_ca(:,:)) 
     3566               zb_sed_ca(:,:) = za_sed_ca(:,:) 
    33953567               zn_sed_ca(:,:) = za_sed_ca(:,:) 
    33963568               za_sed_ca(:,:) = 0.0 
    33973569            endif 
    33983570         ENDIF 
    3399       ELSE                                                   !  case of smolar scheme or muscl 
    3400          if (jorgben.eq.1) then 
    3401             zb_sed_n(:,:)  = za_sed_n(:,:) 
    3402             zn_sed_n(:,:)  = za_sed_n(:,:) 
    3403             za_sed_n(:,:)  = 0.0 
    3404             !! 
    3405             zb_sed_fe(:,:) = za_sed_fe(:,:) 
    3406             zn_sed_fe(:,:) = za_sed_fe(:,:) 
    3407             za_sed_fe(:,:) = 0.0 
    3408             !! 
    3409             zb_sed_c(:,:)  = za_sed_c(:,:) 
    3410             zn_sed_c(:,:)  = za_sed_c(:,:) 
    3411             za_sed_c(:,:)  = 0.0 
    3412          endif 
    3413          if (jinorgben.eq.1) then 
    3414             zb_sed_si(:,:) = za_sed_si(:,:) 
    3415             zn_sed_si(:,:) = za_sed_si(:,:) 
    3416             za_sed_si(:,:) = 0.0 
    3417             !! 
    3418             zb_sed_ca(:,:) = za_sed_ca(:,:) 
    3419             zn_sed_ca(:,:) = za_sed_ca(:,:) 
    3420             za_sed_ca(:,:) = 0.0 
    3421          endif 
    3422       ENDIF 
    34233571      endif 
    3424  
     3572       
    34253573      IF( ln_diatrc ) THEN 
    34263574         !!---------------------------------------------------------------------- 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcco2_medusa.F90

    r5726 r5841  
    4747!======================================================================= 
    4848! 
    49       SUBROUTINE trc_co2_medusa( Temp, Sal, DIC, ALK, Depth, Wnd, pCO2a, & 
     49      SUBROUTINE trc_co2_medusa( Temp, Sal, DIC, ALK, Depth, xkw, pCO2a, & 
    5050      pH, pCO2w, h2co3, hco3, co3, om_cal, om_arg, co2flux, TDIC, TALK,  & 
    5151      dcf, henry, iters ) 
     
    7272! 17/02/2010. Update calculation of K1, K2, Kb to make consistant with the OCMIP protocols. 
    7373! 29/07/2011. Merged into MEDUSA with a raft of changes to this subroutine; less elsewhere 
     74! 23/06/2015. Modified to take gas transfer velocity as an input (rather than wind speed);  
     75!             alter CO2 flux to /s rather than /d for consistency with other schemes 
    7476! 
    7577! Changes for MEDUSA include:  
     
    8587      REAL(wp), INTENT( in )    :: ALK        ! meq  / m3 
    8688      REAL(wp), INTENT( in )    :: Depth      ! m 
    87       REAL(wp), INTENT( in )    :: Wnd        ! m / s 
     89!     REAL(wp), INTENT( in )    :: Wnd        ! m / s 
     90      REAL(wp), INTENT( in )    :: xkw        ! m / s 
    8891      REAL(wp), INTENT( in )    :: pCO2a      ! uatm 
    8992!---------------------------------------------------------------------- 
     
    9598      REAL(wp), INTENT( inout ) :: om_cal     ! normalised 
    9699      REAL(wp), INTENT( inout ) :: om_arg     ! normalised 
    97       REAL(wp), INTENT( inout ) :: co2flux    ! mmol / m2 / d 
     100      REAL(wp), INTENT( inout ) :: co2flux    ! mmol / m2 / s 
    98101      REAL(wp), INTENT( inout ) :: TDIC       ! umol / kg 
    99102      REAL(wp), INTENT( inout ) :: TALK       ! ueq  / kg 
     
    129132   !                 (i.e. surface calculations being performed) 
    130133   if (Depth .eq. 0.0) then 
    131       call Air_sea_exchange( Temp, Wnd, pCO2w, pCO2a, henry, dcf, &         ! inputs 
     134      call Air_sea_exchange( Temp, xkw, pCO2w, pCO2a, henry, dcf, &         ! inputs 
    132135         co2flux )                                                          ! output 
    133136   else 
     
    145148      IF(lwp) WRITE(numout,*) ' trc_co2_medusa: zdic    =', DIC 
    146149      IF(lwp) WRITE(numout,*) ' trc_co2_medusa: zalk    =', ALK 
    147       IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_wind  =', Wnd 
     150      IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_kw660 =', xkw 
    148151      IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_ph    =', ph 
    149152      IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_pco2w =', pCO2w 
     
    162165      &        ' DIC', DIC, ' ALK', ALK 
    163166      if (lwp) write (numout,'(a,a,f10.3,a,f10.3)') 'CO2FLUX-NAN', & 
    164       &        ' WND', Wnd, ' PH ', ph 
     167      &        ' XKW', xkw, ' PH ', ph 
    165168      if (lwp) write (numout,'(a,a,i6)') 'CO2FLUX-NAN', & 
    166169      &        ' ITERS', iters 
     
    196199!  WRITE(*,'(A27,F10.3)') "    Omega calcite       (~) = ", om_cal 
    197200!  WRITE(*,'(A27,F10.3)') "    Omega aragonite     (~) = ", om_arg 
    198 !  WRITE(*,'(A27,F10.3)') "    air sea flux(mmol/m2/d) = ", flux 
     201!  WRITE(*,'(A27,F10.3)') "    air sea flux(mmol/m2/s) = ", flux 
    199202!  WRITE(*,*) " " 
    200203 
     
    287290!======================================================================= 
    288291! 
    289       SUBROUTINE Air_sea_exchange( T, Wnd, pco2w, pco2a, henry, dcf, & 
     292      SUBROUTINE Air_sea_exchange( T, xkw, pco2w, pco2a, henry, dcf, & 
    290293      flux ) 
    291294!       
     
    302305!  pCO2a    partial pressure of CO2 in the atmosphere (usually external forcing). 
    303306!  T        temperature (C) 
    304 !  Wnd      wind speed, metres 
     307!  Wnd      wind speed, metres (DELETED) 
     308!  xkw      gas transfer velocity 
    305309!  Henry    henry's constant 
    306310!  density  the density of water for conversion between mmol/m3 and umol/kg 
     
    312316   IMPLICIT NONE 
    313317 
    314       REAL(wp), INTENT( in )    :: T, wnd, pco2w, pco2a, henry, dcf ! INPUT PARAMETERS: 
     318      REAL(wp), INTENT( in )    :: T, xkw, pco2w, pco2a, henry, dcf ! INPUT PARAMETERS: 
    315319!----------------------------------------------------------------------- 
    316320      REAL(wp), INTENT( inout ) :: flux                             ! OUTPUT Variables 
     
    320324! calculate the Schmidt number and unit conversions 
    321325          sc    = 2073.1-125.62*T+3.6276*T**2.0-0.0432190*T**3.0 
    322           fwind = (0.222d0 * wnd**2d0 + 0.333d0 * wnd)*(sc/660.d0)**(-0.5) 
     326!         fwind = (0.222d0 * wnd**2d0 + 0.333d0 * wnd)*(sc/660.d0)**(-0.5) 
     327          fwind = xkw * (sc/660.d0)**(-0.5) 
    323328          fwind = fwind*24.d0/100.d0   ! convert to m/day 
    324329 
     
    326331! here it is rescaled to mmol/m2/d 
    327332          flux = fwind * henry * ( pco2a - pco2w ) * dcf 
     333 
     334! AXY (23/06/15): let's get it from /d to /s 
     335          flux = flux / ( 86400. ) 
    328336 
    329337  RETURN  
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90

    r5726 r5841  
    3939! 
    4040   SUBROUTINE trc_dms_medusa( chn, chd, mld, xqsr, xdin,   &  !! inputs 
    41      &  dms_surf, dms_andr, dms_simo, dms_aran, dms_hall )    !! outputs 
     41     &  dms_andr, dms_simo, dms_aran, dms_hall )              !! outputs 
    4242!       
    4343!======================================================================= 
     
    7070      !!                    published (and different from the above) 
    7171      !! 
     72      !! AXY (08/07/15): amend to remove Julien's original calculation 
     73      !!                 as this is now superfluous; the four schemes  
     74      !!                 are calculated and one is chosen to be passed 
     75      !!                 to the atmosphere in trc_bio_medusa 
     76      !! 
    7277!======================================================================= 
    7378 
     
    7984      REAL(wp), INTENT( in )    :: xqsr                 !! surface irradiance        (W/m2) 
    8085      REAL(wp), INTENT( in )    :: xdin                 !! surface DIN               (mmol N/m3) 
    81       REAL(wp), INTENT( inout ) :: dms_surf             !! DMS surface concentration (mol/m3)  
    8286      REAL(wp), INTENT( inout ) :: dms_andr             !! DMS surface concentration (mol/m3)  
    8387      REAL(wp), INTENT( inout ) :: dms_simo             !! DMS surface concentration (mol/m3)  
     
    8993      !! temporary variables 
    9094      REAL(wp) ::    fq1,fq2,fq3 
    91 ! 
    92 !! IJT (30/03/13): DMS calc needs this 
    93 !! Julien : in Simo & Dachs, GBC, 2002, DMS is derived from  
    94 !!          CHL/MLD ratio in mg/m4 (i.e. CHL is in mg/m3 
    95 !!          MLD in m). 
    96 !!          In MEDUSA, we already have CHL in mg/m3 for both 
    97 !!          Diatoms and non-diatoms (zchn,zchd); and mld from 
    98 !!          NEMO (hmld) in m. 
    99       CHL = 0.0 
    100 !! 
    101 !!            CHL = mask * TT(I,J,1,PHYTO_TRACER) & 
    102 !!     &       * c2n_p * mw_carbon / CCHL_P(I,J,1,1) 
    103       CHL = chn+chd                                 !! mg/m3  
    104 !! 
    105 !! ------------------------------------------------ 
    106 !!  Calculate the DMS concentration in nM (nanomol/litre) 
    107 !!   from Simo & Dachs, GBC, 2002, modified to be positive-definite 
    108 !!   for MLD>182.536m, using DMS=90./MLD (Aranami & Tsunogai, JGR, 2004) 
    109 !!   Multiply by 1.0E-6 to convert nM to (mol/m3) 
    110 !!       cmr = fm(i,1)*chl/mld(i) 
    111 !!       IF (cmr .lt. 0.02) THEN 
    112 !!         IF (mld(i) .le. 182.536) THEN 
    113 !!           csurf(i,astf_dms) = 1.0e-6*fm(i,1)*(-LOG(mld(i)) + 5.7) 
    114 !!         ELSE 
    115 !!           csurf(i,astf_dms) = 1.0e-6*fm(i,1)*(90./mld(i)) 
    116 !!         ENDIF 
    117 !!       ELSE 
    118 !!         csurf(i,astf_dms) = 1.0e-6*fm(i,1)*(55.8*cmr + 0.6) 
    119 !!       ENDIF 
    120 !! 
    121         cmr      = CHL / mld 
    122 !       sw_dms   = 0.5 + SIGN( 0.5, cmr - 0.02 ) 
    123 !! Jpalm (11-08-2014) 
    124 !! Explanation about the SIGN function : 
    125 !! not easy to read, but maybe "more elegant and efficient") 
    126 !! here for example:  
    127 !! sw_dms = 1 if cmr is greater than 0.02, 
    128 !!          0 if cmr lower than 0.02 
    129 !! then  
    130 !! if cmr < 0.02 
    131 !!  dms_surf =  1.0e-6 * 90.0 / mld  
    132 !!       or  =  1.0e-6 * 5.7 - LOG(mld) 
    133 !! and if cmr > 0.02 
    134 !!  dms_surf = 1.0e-6 * ( 55.8 * cmr + 0.6 ) 
    135 !! what is equivalent to the IF loops formulations. 
    136 !! difference is on the stresholds between mld = 182.536m 
    137 !! (strange value...) 
    138 !! and the Max function... that stay uncertain. 
    139 !! 
    140 !        dms_surf = 1.0e-6 * ( sw_dms *             & 
    141 !     &  ( 55.8 * cmr + 0.6 ) + ( 1.0 - sw_dms ) *  & 
    142 !     &  ( MAX( 90.0 / mld, 5.7 - LOG(mld) ) ) ) 
    143 ! 
    144 ! AXY (12/01/15): the DMS equation donated by the UKMO does not match 
    145 !                 that reported in Halloran et al. (2010); amend the 
    146 !                 equations appropriately 
    147 ! 
    148         if (cmr .lt. 0.02) then 
    149            dms_surf = (-1.0 * log(mld)) + 5.7 
    150         else 
    151            dms_surf = (55.8 * cmr) + 0.6 
    152         endif 
    153 !     
    154         if (mld > 182.5) then 
    155            dms_surf = (90.0 / mld) 
    156         endif 
    157 !      
    158         dms_surf = 1.0e-6 * dms_surf 
    159  
    16095!  
    16196!======================================================================= 
     
    16398! AXY (13/03/15): per remarks above, the following calculations estimate 
    16499!                 DMS using all of the schemes examined for UKESM1 
     100! 
     101      CHL = 0.0 
     102      CHL = chn+chd                                 !! mg/m3  
     103      cmr = CHL / mld 
    165104! 
    166105! AXY (13/03/15): Anderson et al. (2001) 
     
    180119! 
    181120! AXY (13/03/15): Simo & Dachs (2002) 
    182         cmr = CHL / mld 
    183121        fq1 = (-1 * log(mld)) + 5.7 
    184122        fq2 = (55.8 * cmr) + 0.6 
     
    191129!            
    192130! AXY (13/03/15): Aranami & Tsunogai (2004) 
    193         cmr = CHL / mld 
    194131        fq1 = 60.0 / mld 
    195132        fq2 = (55.8 * cmr) + 0.6 
     
    202139!         
    203140! AXY (13/03/15): Halloran et al. (2010) 
    204         cmr = CHL / mld 
    205141        fq1 = (-1 * log(mld)) + 5.7 
    206142        fq2 = (55.8 * cmr) + 0.6 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90

    r5741 r5841  
    3838   !! JPALM (14/09/15) 
    3939   LOGICAL, PUBLIC ::                  & 
    40       ln_ccd = .FALSE. 
     40      ln_ccd = .TRUE. 
    4141 
    4242   INTEGER ::                          & 
     
    236236      
    237237      !!---------------------------------------------------------------------- 
     238      !! Averaged properties for DMS calculations (various units) 
     239      !!---------------------------------------------------------------------- 
     240      !!      
     241      !! these store temporally averaged properties for DMS calculations (AXY, 07/07/15) 
     242      zb_dms_chn(:,:)  = 0.0  !! CHN 
     243      zn_dms_chn(:,:)  = 0.0 
     244      za_dms_chn(:,:)  = 0.0 
     245      zb_dms_chd(:,:)  = 0.0  !! CHD 
     246      zn_dms_chd(:,:)  = 0.0 
     247      za_dms_chd(:,:)  = 0.0 
     248      zb_dms_mld(:,:)  = 0.0  !! MLD 
     249      zn_dms_mld(:,:)  = 0.0 
     250      za_dms_mld(:,:)  = 0.0 
     251      zb_dms_qsr(:,:)  = 0.0  !! QSR 
     252      zn_dms_qsr(:,:)  = 0.0 
     253      za_dms_qsr(:,:)  = 0.0 
     254      zb_dms_din(:,:)  = 0.0  !! DIN 
     255      zn_dms_din(:,:)  = 0.0 
     256      za_dms_din(:,:)  = 0.0 
     257      !! 
     258      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: average fields for DMS initialised to zero' 
     259      IF(lwp) CALL flush(numout) 
     260 
     261      !!---------------------------------------------------------------------- 
    238262      !! AXY (04/11/13): initialise fields previously done by trc_sed_medusa 
    239263      !!---------------------------------------------------------------------- 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90

    r5736 r5841  
    8383      &  xsdiss,                                              & 
    8484      &  vsed,xhr,                                            & 
    85       &  sedlam,sedlostpoc,jpkb,jdms 
     85      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model 
    8686#if defined key_roam 
    8787      NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        & 
     
    138138     IF( ( .NOT.lk_iomput .AND. ln_diatrc ) .OR. ( ln_diatrc .AND. lk_medusa ) ) THEN 
    139139         ! 
    140          ! Namelist nampisdia 
     140         ! Namelist nammeddia 
    141141         ! ------------------- 
    142          REWIND( numnatp_ref )              ! Namelist nampisdia in reference namelist : Pisces diagnostics 
     142         REWIND( numnatp_ref )              ! Namelist nammeddia in reference namelist : MEDUSA diagnostics 
    143143         READ  ( numnatp_ref, nammeddia, IOSTAT = ios, ERR = 901) 
    144144901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in reference namelist', lwp ) 
    145145 
    146          REWIND( numnatp_cfg )              ! Namelist nampisdia in configuration namelist : Pisces diagnostics 
     146         REWIND( numnatp_cfg )              ! Namelist nammeddia in configuration namelist : MEDUSA diagnostics 
    147147         READ  ( numnatp_cfg, nammeddia, IOSTAT = ios, ERR = 902 ) 
    148148902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in configuration namelist', lwp ) 
     
    338338      jpkb    = 0. 
    339339      jdms        = 0 
     340      jdms_input  = 0 
     341      jdms_input  = 3 
    340342             
    341343      !REWIND(numnatm) 
     
    343345         ! Namelist natbio 
    344346         ! ------------------- 
    345          REWIND( numnatp_ref )              ! Namelist nampisdia in reference namelist : Pisces diagnostics 
     347         REWIND( numnatp_ref )              ! Namelist natbio in reference namelist : MEDUSA diagnostics 
    346348         READ  ( numnatp_ref, natbio, IOSTAT = ios, ERR = 903) 
    347 903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in reference namelist', lwp ) 
    348  
    349          REWIND( numnatp_cfg )              ! Namelist nampisdia in configuration namelist : Pisces diagnostics 
     349903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natbio in reference namelist', lwp ) 
     350 
     351         REWIND( numnatp_cfg )              ! Namelist natbio in configuration namelist : MEDUSA diagnostics 
    350352         READ  ( numnatp_cfg, natbio, IOSTAT = ios, ERR = 904 ) 
    351 904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in configuration namelist', lwp ) 
     353904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natbio in configuration namelist', lwp ) 
    352354         IF(lwm) WRITE ( numonp, natbio ) 
    353355 
     
    488490!! UKESM1 - new diagnostics  !! Jpalm 
    489491!!       jdms        :  include dms diagnostics 
    490 !! 
    491 !! 
    492 !! 
    493  
     492!!  jdms_input  :  use instant (0) or diel-avg (1) inputs 
     493!!       jdms_model  :  choice of DMS model passed to atmosphere 
     494!!                      1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL 
     495!! 
    494496      IF(lwp) THEN 
    495497!! 
     
    510512         &   ' key_axy_carbchem                                                       = INACTIVE' 
    511513#endif         
     514#if defined key_mocsy 
     515         WRITE(numout,*)     & 
     516         &   ' key_mocsy                                                              = ACTIVE' 
     517#else 
     518         WRITE(numout,*)     & 
     519         &   ' key_mocsy                                                              = INACTIVE' 
     520#endif         
     521#if defined key_avgqsr_medusa 
     522         WRITE(numout,*)     & 
     523         &   ' key_avgqsr_medusa                                                      = ACTIVE' 
     524#else 
     525         WRITE(numout,*)     & 
     526         &   ' key_avgqsr_medusa                                                      = INACTIVE' 
     527#endif         
    512528#if defined key_bs_axy_zforce 
    513529         WRITE(numout,*)     & 
     
    544560         WRITE(numout,*)     & 
    545561         &   ' key_axy_pi_co2                                                         = INACTIVE' 
     562# endif 
     563# if defined key_debug_medusa 
     564         WRITE(numout,*)     & 
     565         &   ' key_debug_medusa                                                       = ACTIVE' 
     566#else 
     567         WRITE(numout,*)     & 
     568         &   ' key_debug_medusa                                                       = INACTIVE' 
    546569# endif 
    547570         WRITE(numout,*) ' ' 
     
    971994         &   ' Vert layer for diagnostic of vertical flux,                jpkp        = ', jpkb 
    972995!! 
    973 !! UKESM1 - new diagnostics  !! Jpalm 
     996!! UKESM1 - new diagnostics  !! Jpalm; AXY (08/07/15) 
    974997         WRITE(numout,*) '=== UKESM1-related parameters' 
    975998         WRITE(numout,*)     & 
    976999         &   ' include DMS diagnostic?,                                   jdms        = ', jdms 
     1000         if (jdms_input .eq. 0) then 
     1001            WRITE(numout,*)     & 
     1002            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous' 
     1003         else 
     1004            WRITE(numout,*)     & 
     1005            &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average' 
     1006         endif 
     1007    if (jdms_model .eq. 1) then 
     1008            WRITE(numout,*)     & 
     1009            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)' 
     1010    elseif (jdms_model .eq. 2) then 
     1011            WRITE(numout,*)     & 
     1012            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)' 
     1013    elseif (jdms_model .eq. 3) then 
     1014            WRITE(numout,*)     & 
     1015            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)' 
     1016    elseif (jdms_model .eq. 4) then 
     1017            WRITE(numout,*)     & 
     1018            &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)' 
     1019         endif 
    9771020!! 
    9781021      ENDIF 
     
    10321075 
    10331076      !READ(numnatm,natroam) 
    1034          ! Namelist natbio 
     1077         ! Namelist natroam 
    10351078         ! ------------------- 
    1036          REWIND( numnatp_ref )              ! Namelist nampisdia in reference namelist : Pisces diagnostics 
    1037          READ  ( numnatp_ref, natbio, IOSTAT = ios, ERR = 905) 
    1038 905      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in reference namelist', lwp ) 
    1039  
    1040          REWIND( numnatp_cfg )              ! Namelist nampisdia in configuration namelist : Pisces diagnostics 
    1041          READ  ( numnatp_cfg, natbio, IOSTAT = ios, ERR = 906 ) 
    1042 906      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in configuration namelist', lwp ) 
    1043          IF(lwm) WRITE ( numonp, natbio ) 
     1079         REWIND( numnatp_ref )              ! Namelist natroam in reference namelist : MEDUSA diagnostics 
     1080         READ  ( numnatp_ref, natroam, IOSTAT = ios, ERR = 905) 
     1081905      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natroam in reference namelist', lwp ) 
     1082 
     1083         REWIND( numnatp_cfg )              ! Namelist natroam in configuration namelist : MEDUSA diagnostics 
     1084         READ  ( numnatp_cfg, natroam, IOSTAT = ios, ERR = 906 ) 
     1085906      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natroam in configuration namelist', lwp ) 
     1086         IF(lwm) WRITE ( numonp, natroam ) 
    10441087 
    10451088!! ROAM carbon, alkalinity and oxygen cycle parameters 
     
    10861129         ! Namelist natopt 
    10871130         ! ------------------- 
    1088          REWIND( numnatp_ref )              ! Namelist nampisdia in reference namelist : Pisces diagnostics 
     1131         REWIND( numnatp_ref )              ! Namelist natopt in reference namelist : MEDUSA diagnostics 
    10891132         READ  ( numnatp_ref, natopt, IOSTAT = ios, ERR = 907) 
    1090 907      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in reference namelist', lwp ) 
    1091  
    1092          REWIND( numnatp_cfg )              ! Namelist nampisdia in configuration namelist : Pisces diagnostics 
     1133907      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natopt in reference namelist', lwp ) 
     1134 
     1135         REWIND( numnatp_cfg )              ! Namelist natopt in configuration namelist : MEDUSA diagnostics 
    10931136         READ  ( numnatp_cfg, natopt, IOSTAT = ios, ERR = 908 ) 
    1094 908      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammeddia in configuration namelist', lwp ) 
     1137908      IF( ios /= 0 ) CALL ctl_nam ( ios , 'natopt in configuration namelist', lwp ) 
    10951138         IF(lwm) WRITE ( numonp, natopt ) 
    10961139 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90

    r5726 r5841  
    7373      ! determination of surface irradiance 
    7474      ! ----------------------------------- 
    75       zpar0m (:,:)   = qsr   (:,:) * 0.43 
     75      ! AXY (23/07/15): the inclusion of empirical DMS calculations requires 
     76      !                 daily averages of a series of properties that are 
     77      !                 used as inputs; these include surface irradiance;  
     78      !                 here, this is taken advantage of to allow MEDUSA to 
     79      !                 base its submarine light field on daily average 
     80      !                 rather than "instantaneous" irradiance; largely 
     81      !                 because MEDUSA was originally formulated to work 
     82      !                 with diel average irradiance rather than a diel 
     83      !                 cycle; using key_avgqsr_medusa activates this 
     84      !                 functionality, while its absence gives default 
     85      !                 MEDUSA (which is whatever is supplied by NEMO) 
     86# if defined key_avgqsr_medusa 
     87      ! surface irradiance input is rolling average irradiance 
     88      zpar0m (:,:)   = zn_dms_qsr(:,:) * 0.43 
     89# else       
     90      ! surface irradiance input is   instantaneous irradiance 
     91      zpar0m (:,:)   =        qsr(:,:) * 0.43 
     92# endif 
    7693      ! AXY (22/08/14): when zpar0m = 0, zpar100 is also zero and calculating  
    7794      !                 euphotic depth is not possible (cf. the Arctic Octopus);  
     
    92109      zparr  (:,:,1) = 0.5 * zpar0m(:,:) 
    93110      zparg  (:,:,1) = 0.5 * zpar0m(:,:) 
    94  
    95111 
    96112      ! determination of xpar 
     
    146162      ENDDO  
    147163 
    148  
    149164      IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 
    150165         WRITE(charout, FMT="('opt')") 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcoxy_medusa.F90

    r5726 r5841  
    3434! The following is a map of the subroutines contained within this module 
    3535! - trc_oxy_medusa 
    36 !      - CALLS gas_transfer 
    3736!      - CALLS oxy_schmidt 
    3837!      - CALLS oxy_sato 
     
    4241!======================================================================= 
    4342! 
    44    SUBROUTINE trc_oxy_medusa( pt, ps, uwind, vwind, pp0, o2, dz, &  !! inputs 
    45       kw660, o2flux, o2sat )                                        !! outputs 
     43   SUBROUTINE trc_oxy_medusa( pt, ps, kw660, pp0, o2, &  !! inputs 
     44      kwo2, o2flux, o2sat )                               !! outputs 
    4645!       
    4746!======================================================================= 
     
    5756      !! number of oxygen) and oxy_sato.f (calculates oxygen saturation  
    5857      !! concentration at 1 atm). 
     58      !! 
     59      !! AXY (23/06/15): revised to allow common gas transfer velocity 
     60      !!                 to be used for CO2 and O2; outputs of this 
     61      !!                 routine amended to mmol/m3 from mol/m3 
    5962      !!  
    6063      !! Function inputs are (in order) :  
    6164      !!     pt      temperature                     (degrees C) 
    6265      !!     ps      salinity                        (PSU) 
    63       !!     uwind   u-wind velocity                 (m/s) 
    64       !!     vwind   v-wind velocity                 (m/s) 
     66      !!     kw660   gas transfer velocity           (m/s) 
    6567      !!     pp0     surface pressure                (divided by 1 atm) 
    66       !!     o2      surface O2 concentration        (mol/m3) 
    67       !!     dz      surface layer thickness         (m) 
    68       !! (*) kw660   gas transfer velocity           (m/s) 
    69       !! (*) o2flux  exchange rate of oxygen         (mol/m3/s) 
    70       !! (+) o2sat   oxygen saturation concentration (mol/m3) 
     68      !!     o2      surface O2 concentration        (mmol/m3) 
     69      !! (+) kwo2    gas transfer velocity for O2    (m/s) 
     70      !! (*) o2flux  exchange rate of oxygen         (mmol/m2/s) 
     71      !! (+) o2sat   oxygen saturation concentration (mmol/m3) 
    7172      !!  
    7273      !! Where (*) is the function output (note its units).   
     
    7879      REAL(wp), INTENT( in )    :: pt 
    7980      REAL(wp), INTENT( in )    :: ps 
    80       REAL(wp), INTENT( in )    :: uwind 
    81       REAL(wp), INTENT( in )    :: vwind 
     81      REAL(wp), INTENT( in )    :: kw660 
    8282      REAL(wp), INTENT( in )    :: pp0 
    8383      REAL(wp), INTENT( in )    :: o2 
    84       REAL(wp), INTENT( in )    :: dz 
    85       REAL(wp), INTENT( inout ) :: kw660, o2flux, o2sat 
    86 ! 
    87       REAL(wp) :: scl_wind, kwo2, o2schmidt, o2sato 
    88 ! 
    89 ! Calculate gas transfer 
    90 !  
    91       call gas_transfer(uwind, vwind, scl_wind, kw660) 
     84      REAL(wp), INTENT( out )   :: kwo2, o2flux, o2sat 
     85! 
     86      REAL(wp) :: o2schmidt, o2sato, mol_o2 
     87! 
     88! Oxygen to mol / m3 
     89! 
     90      mol_o2 = o2 / 1000. 
    9291! 
    9392! Calculate oxygen Schmidt number 
     
    106105! Calculate time rate of change of O2 due to gas exchange (mol/m3/s) 
    107106! 
    108       o2flux = kwo2 * (o2sat - o2) / dz 
     107      o2flux = kwo2 * (o2sat - mol_o2) 
     108! 
     109! Oxygen flux and saturation to mmol / m3 
     110! 
     111      o2sat  =  o2sat * 1000. 
     112      o2flux = o2flux * 1000. 
    109113! 
    110114      END SUBROUTINE trc_oxy_medusa 
     
    130134      !! are taken from Keeling et al. (1998, GBC, 12, 141-163). 
    131135      !!  
     136      !! AXY (23/06/2015) 
     137      !! UPDATED: revised formulation from Wanninkhof (2014) for 
     138      !! consistency with MOCSY 
     139      !! 
     140      !! Winninkhof, R. (2014). Relationship between wind speed and gas 
     141      !! exchange over the ocean revisited. LIMNOLOGY AND OCEANOGRAPHY-METHODS 
     142      !! 12, 351-362, doi:10.4319/lom.2014.12.351 
     143      !! 
    132144      !! Function inputs are (in order) :  
    133145      !!     t           temperature (degrees C) 
     
    141153! 
    142154      REAL(wp) :: pt, o2_schmidt 
    143       REAL(wp) :: a0, a1, a2, a3 
    144 ! 
    145       data a0 /    1638.0 / 
    146       data a1 /    -81.83 / 
    147       data a2 /     1.483 / 
    148       data a3 / -0.008004 / 
    149 ! 
    150       o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*a3)) 
     155      REAL(wp) :: a0, a1, a2, a3, a4 
     156! 
     157! AXY (23/06/15): OCMIP-2 coefficients 
     158!     data a0 /    1638.0 / 
     159!     data a1 /    -81.83 / 
     160!     data a2 /     1.483 / 
     161!     data a3 / -0.008004 / 
     162! 
     163! AXY (23/06/15): Wanninkhof (2014) coefficients 
     164      data a0 /     1920.4 / 
     165      data a1 /     -135.6 / 
     166      data a2 /     5.2121 / 
     167      data a3 /   -0.10939 / 
     168      data a4 / 0.00093777 / 
     169! 
     170!     o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*a3)) 
     171      o2_schmidt = a0 + pt*(a1 + pt*(a2 + pt*(a3 + pt*a4))) 
    151172! 
    152173      END SUBROUTINE oxy_schmidt 
     
    231252!======================================================================= 
    232253 
    233 !======================================================================= 
    234 ! 
    235    SUBROUTINE gas_transfer( uwind, vwind, &  !! input  
    236       scl_wind, k )                          !! output 
    237 !       
    238 !======================================================================= 
    239       !! 
    240       !! Title  : Calculates gas transfer velocity 
    241       !! Author : Andrew Yool 
    242       !! Date   : 15/10/04 (revised 04/08/2011) 
    243       !!  
    244       !! This subroutine uses near-surface wind speed to calculate gas 
    245       !! transfer velocity for use in CO2 and O2 exchange calculations. 
    246       !!  
    247       !! Note that the parameterisation of Wanninkhof quoted here is a 
    248       !! truncation of the original equation.  It excludes a chemical 
    249       !! enhancement function (based on temperature), although such 
    250       !! temperature dependence is reported negligible by Etcheto & 
    251       !! Merlivat (1988). 
    252       !!  
    253       !! Note also that in calculating scalar wind, the variance of the 
    254       !! wind over the period of a timestep is ignored.  Some authors, 
    255       !! for instance OCMIP-2, favour including some reference to the 
    256       !! variability of wind.  However, their wind fields are averaged 
    257       !! over relatively long time periods, and so this issue may be 
    258       !! safely (!) ignored here. 
    259       !!  
    260       !! Subroutine inputs are (in order) : 
    261       !!     uwind     wind u velocity at 10 m (m/s) 
    262       !!     vwind     wind v velocity at 10 m (m/s) 
    263       !! (+) scl_wind  scalar wind velocity at 10 m (m/s) 
    264       !! (*) k         gas transfer velocity (m/s) 
    265       !! Where (*) is the function output and (+) is a diagnostic output. 
    266       !! 
    267 !======================================================================= 
    268  
    269       implicit none 
    270 ! 
    271 ! Input variables 
    272       REAL(wp) :: uwind, vwind 
    273 ! 
    274 ! Output variables 
    275       REAL(wp) :: scl_wind, k, tmp_k 
    276 ! 
    277 ! Choice of parameterisation 
    278       INTEGER  :: eqn 
    279 ! 
    280 ! Coefficients for various parameterisations 
    281       REAL(wp), DIMENSION(6) :: a 
    282       REAL(wp), DIMENSION(6) :: b 
    283 ! 
    284 ! Values of coefficients 
    285       data a(1) / 0.166 /  ! Liss & Merlivat (1986)    [approximated] 
    286       data a(2) / 0.3   /  ! Wanninkhof (1992)         [sans enhancement] 
    287       data a(3) / 0.23  /  ! Nightingale et al. (2000) [good] 
    288       data a(4) / 0.23  /  ! Nightingale et al. (2000) [better] 
    289       data a(5) / 0.222 /  ! Nightingale et al. (2000) [best] 
    290       data a(6) / 0.337 /  ! OCMIP-2                   [sans variability] 
    291 ! 
    292       data b(1) / 0.133 / 
    293       data b(2) / 0.0   / 
    294       data b(3) / 0.0   / 
    295       data b(4) / 0.1   / 
    296       data b(5) / 0.333 / 
    297       data b(6) / 0.0   / 
    298 ! 
    299 ! Which parameterisation is to be used? 
    300       eqn = 2 
    301 ! 
    302 ! Calculate scalar wind (m/s) 
    303       scl_wind = (uwind**2 + vwind**2)**0.5 
    304 ! 
    305 ! Calculate gas transfer velocity (cm/h) 
    306       tmp_k = (a(eqn) * scl_wind**2) + (b(eqn) * scl_wind) 
    307 ! 
    308 ! Convert tmp_k from cm/h to m/s 
    309       k = tmp_k / (100. * 3600.) 
    310 ! 
    311       END SUBROUTINE gas_transfer 
    312  
    313 !======================================================================= 
    314 !======================================================================= 
    315 !======================================================================= 
    316  
    317254#else 
    318255   !!====================================================================== 
     
    322259CONTAINS 
    323260 
    324    SUBROUTINE trc_oxy_medusa( pt, ps, uwind, vwind, pp0, o2, dz, &  !! inputs 
    325       kw660, o2flux, o2sat )                                        !! outputs 
     261   SUBROUTINE trc_oxy_medusa( pt, ps, kw660, pp0, o2, &  !! inputs 
     262      o2flux, o2sat )                                     !! outputs 
    326263      USE par_kind 
    327264 
    328265      REAL(wp), INTENT( in )    :: pt 
    329266      REAL(wp), INTENT( in )    :: ps 
     267      REAL(wp), INTENT( in )    :: kw660 
    330268      REAL(wp), INTENT( in )    :: pp0 
    331269      REAL(wp), INTENT( in )    :: o2 
    332       REAL(wp), INTENT( in )    :: dz 
    333       REAL(wp), INTENT( inout ) :: kw660, o2flux, o2sat 
     270      REAL(wp), INTENT( inout ) :: o2flux, o2sat 
    334271 
    335272      WRITE(*,*) 'trc_oxy_medusa: You should not have seen this print! error?', kt 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsed_medusa.F90

    r5742 r5841  
    4040   !! AXY (10/02/09) 
    4141   LOGICAL, PUBLIC ::                  & 
    42       bdustfer = .FALSE. 
     42      bdustfer = .TRUE. 
    4343   REAL(wp), PUBLIC ::                 & 
    4444      sedfeinput = 1.e-9_wp  ,         & 
     
    322322               CALL prihre( dustmo(:,:,1),jpi,jpj,1,jpi,20,1,jpj,10,1e9,numout ) 
    323323            ENDIF 
    324   
     324 
    325325         ENDIF 
    326326 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90

    r5726 r5841  
    2121   USE trcopt_medusa 
    2222   USE trcsed_medusa 
     23   USE trcavg_medusa 
    2324 
    2425 
     
    4647      INTEGER, INTENT(in) :: kt   ! ocean time-step index 
    4748 
     49# if defined key_debug_medusa 
     50         IF(lwp) WRITE(numout,*) ' MEDUSA inside trc_sms_medusa' 
     51         CALL flush(numout) 
     52# endif 
     53 
    4854      IF( kt == nittrc000 ) THEN 
    4955       IF(lwp) WRITE(numout,*) 
     
    5258      ENDIF 
    5359 
     60      CALL trc_avg_medusa( kt )   ! rolling average module 
     61# if defined key_debug_medusa 
     62         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_avg_medusa' 
     63         CALL flush(numout) 
     64# endif 
     65 
    5466      CALL trc_opt_medusa( kt )   ! optical model 
     67# if defined key_debug_medusa 
     68         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_opt_medusa' 
     69         CALL flush(numout) 
     70# endif 
    5571 
    5672# if defined key_kill_medusa 
     
    6076# else 
    6177      CALL trc_bio_medusa( kt )   ! biological model 
     78# if defined key_debug_medusa 
     79         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_bio_medusa' 
     80         CALL flush(numout) 
     81# endif 
    6282 
    6383      CALL trc_sed_medusa( kt )   ! sedimentation model 
     84# if defined key_debug_medusa 
     85         IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa' 
     86         CALL flush(numout) 
     87# endif 
    6488# endif 
    6589 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/par_trc.F90

    r5726 r5841  
    88   !!             1.0  !  2004-03  (C. Ethe) Free form and module 
    99   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  revised architecture 
     10   !!              -   !  2014-06  (A. Yool, J. Palmieri) adding MEDUSA-2 
    1011   !!---------------------------------------------------------------------- 
    1112   USE par_kind          ! kind parameters 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcnam.F90

    r5726 r5841  
    1111   !!              -   !  2001-01 (E Kestenare) suppress ndttrc=1 for CEN2 and TVD schemes 
    1212   !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90 
     13   !!              -   !  2014-06  (A. Yool, J. Palmieri) adding MEDUSA-2 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_top 
     
    5657      !! ** Method  : - read passive tracer namelist  
    5758      !!              - read namelist of each defined SMS model 
    58       !!                ( (PISCES, CFC, MY_TRC ) 
     59      !!                ( (PISCES, CFC, MY_TRC, MEDUSA, IDTRA ) 
    5960      !!--------------------------------------------------------------------- 
    6061      INTEGER  ::   jn, jk                     ! dummy loop indice 
     
    231232      ENDIF 
    232233      ! 
    233 # if defined key_debug_medusa 
    234       CALL flush(numout) 
    235       IF (lwp) write (numout,*) '------------------------------' 
    236       IF (lwp) write (numout,*) 'Jpalm - debug' 
    237       IF (lwp) write (numout,*) 'CALL trc_nam_pisces -- OK' 
     234 
     235      IF( lk_cfc     ) THEN   ;   CALL trc_nam_cfc         ! CFC     tracers 
     236      ELSE                    ;   IF(lwp) WRITE(numout,*) '          CFC not used' 
     237      ENDIF 
     238 
     239      IF( lk_c14b     ) THEN   ;   CALL trc_nam_c14b         ! C14 bomb     tracers 
     240      ELSE                    ;   IF(lwp) WRITE(numout,*) '          C14 not used' 
     241      ENDIF 
     242 
     243      IF( lk_my_trc  ) THEN   ;   CALL trc_nam_my_trc      ! MY_TRC  tracers 
     244      ELSE                    ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used' 
     245      ENDIF 
     246      ! 
     247# if defined key_debug_medusa 
     248      CALL flush(numout) 
     249      IF (lwp) write (numout,*) '------------------------------' 
     250      IF (lwp) write (numout,*) 'Jpalm - debug' 
     251      IF (lwp) write (numout,*) 'CALL trc_nam_pisces - CFC - C14 - my_trc -- OK' 
    238252      IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_medusa' 
    239253      IF (lwp) write (numout,*) ' ' 
     
    262276      IF (lwp) write (numout,*) 'Jpalm - debug' 
    263277      IF (lwp) write (numout,*) 'CALL trc_nam_idtra -- OK' 
    264       IF (lwp) write (numout,*) 'in trc_nam - just before CALL trc_nam_cfc' 
    265       IF (lwp) write (numout,*) ' ' 
    266 # endif 
    267       ! 
    268       IF( lk_cfc     ) THEN   ;   CALL trc_nam_cfc         ! CFC     tracers 
    269       ELSE                    ;   IF(lwp) WRITE(numout,*) '          CFC not used' 
    270       ENDIF 
    271  
    272       IF( lk_c14b     ) THEN   ;   CALL trc_nam_c14b         ! C14 bomb     tracers 
    273       ELSE                    ;   IF(lwp) WRITE(numout,*) '          C14 not used' 
    274       ENDIF 
    275  
    276       IF( lk_my_trc  ) THEN   ;   CALL trc_nam_my_trc      ! MY_TRC  tracers 
    277       ELSE                    ;   IF(lwp) WRITE(numout,*) '          MY_TRC not used' 
    278       ENDIF 
     278      IF (lwp) write (numout,*) 'in trc_nam - CALL trc_nam OK' 
     279      IF (lwp) write (numout,*) ' ' 
     280# endif 
    279281      ! 
    280282      IF(lwp)   CALL flush(numout) 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcsms.F90

    r5726 r5841  
    5151      ! 
    5252      IF( lk_pisces  )   CALL trc_sms_pisces ( kt )    ! main program of PISCES  
    53       IF( lk_medusa  )   CALL trc_sms_medusa ( kt )    ! MEDUSA  tracers 
    54       IF( lk_idtra   )   CALL trc_sms_idtra  ( kt )    ! radioactive decay of Id. tracer 
    5553      IF( lk_cfc     )   CALL trc_sms_cfc    ( kt )    ! surface fluxes of CFC 
    5654      IF( lk_c14b    )   CALL trc_sms_c14b   ( kt )    ! surface fluxes of C14 
    5755      IF( lk_my_trc  )   CALL trc_sms_my_trc ( kt )    ! MY_TRC  tracers 
     56      IF( lk_medusa  )   CALL trc_sms_medusa ( kt )    ! MEDUSA  tracers 
     57      IF( lk_idtra   )   CALL trc_sms_idtra  ( kt )    ! radioactive decay of Id. tracer 
    5858 
    5959      IF(ln_ctl) THEN      ! print mean trends (used for debugging) 
  • branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/trcwri.F90

    r5726 r5841  
    55   !!====================================================================== 
    66   !! History :   1.0  !  2009-05 (C. Ethe)  Original code 
     7   !!              -   !  2014-06 (A. Yool, J. Palmieri) adding MEDUSA-2 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_top && defined key_iomput 
     
    5859      ! --------------------------------------- 
    5960      IF( lk_pisces  )   CALL trc_wri_pisces     ! PISCES  
     61      IF( lk_cfc     )   CALL trc_wri_cfc        ! surface fluxes of CFC 
     62      IF( lk_c14b    )   CALL trc_wri_c14b       ! surface fluxes of C14 
     63      IF( lk_my_trc  )   CALL trc_wri_my_trc     ! MY_TRC  tracers 
    6064      ! 
    6165# if defined key_debug_medusa 
     
    8185      !!! JPALM 
    8286      !!! don't forget to add idtra  
    83       IF( lk_cfc     )   CALL trc_wri_cfc        ! surface fluxes of CFC 
    84       IF( lk_c14b    )   CALL trc_wri_c14b       ! surface fluxes of C14 
    85       IF( lk_my_trc  )   CALL trc_wri_my_trc     ! MY_TRC  tracers 
    8687      ! 
    8788      IF( nn_timing == 1 )  CALL timing_stop('trc_wri') 
Note: See TracChangeset for help on using the changeset viewer.