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 for branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 – NEMO

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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         !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.