Changeset 8012


Ignore:
Timestamp:
2017-05-10T09:47:35+02:00 (3 years ago)
Author:
marc
Message:

Pulled further code from trcbio_medusa.F90 into other files

Location:
branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90

    r7998 r8012  
    5050                                   fccd, fcomm_resp,                     & 
    5151                                   fd_cal3, fd_car3, fd_nit3, fd_sil3,   & 
     52                                   fdep1,                                & 
    5253                                   fdd, fdd2d, fddc,                     & 
    5354                                   fdpd, fdpd2, fdpd22d, fdpd2d,         & 
     
    7475                                   fprds, fprds2d,                       &  
    7576                                   fprn, fprn_ml, fprn2d,                & 
    76                                    fregen2d,                     & 
     77                                   fregen, fregen2d,                     & 
    7778                                   fregenfast, fregenfastsi,             & 
    78                                    fregensi2d,                 & 
     79                                   fregensi, fregensi2d,                 & 
    7980                                   freminc, freminc2d,                   & 
    8081                                   freminca, freminca2d,                 &  
     
    9596                                   iben_c2d, iben_ca2d, iben_fe2d,       & 
    9697                                   iben_n2d, iben_si2d,                  & 
     98                                   intdissic, intdissin,                 & 
     99                                   intdissisi, inttalk,                  & 
    97100                                   iters, lyso_ca2d,                     & 
    98101                                   mdetc2d, megrazd3, megrazp3,          & 
    99102                                   megrazz3,                             &  
    100103                                   migrazd3, migrazp3,                   & 
     104                                   o2min, o2sat3,                        & 
    101105                                   oben_c2d, oben_ca2d, oben_fe2d,       & 
    102106                                   oben_n2d, oben_si2d,                  & 
     
    114118                                   zimesd2d, zimesdc2d, zimesn2d,        & 
    115119                                   ziresp2d,                             & 
    116                                    zpds, zphd, zphn 
    117       USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask 
     120                                   zo2min,                               & 
     121                                   zalk, zdet, zdic, zdin, zdtc,         & 
     122                                   zfer, zoxy, zpds, zphd, zphn,         & 
     123                                   zsal, zsil, ztmp, zzme, zzmi 
     124      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,       & 
     125                                   mbathy, tmask 
    118126      USE in_out_manager,    ONLY: lwp, numout 
    119127      USE iom,               ONLY: lk_iomput 
    120128      USE par_kind,          ONLY: wp 
    121       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1 
     129      USE par_oce,           ONLY: jpim1, jpjm1 
    122130      USE phycst,            ONLY: rsmall 
    123131      USE sbc_oce,           ONLY: qsr, wndm 
     
    126134                                   i0100, i0150, i0200, i0500, i1000,    & 
    127135                                   jdms, ocal_ccd,                       & 
    128                                    xbetac, xbetan, xpar, xphi,           & 
    129                                    xthetapd, xthetapn, xthetazmi, xze 
     136                                   xbetac, xbetan, xpar, xphi, xrfn,     & 
     137                                   xthetapd, xthetapn, xthetazme, xthetazmi, xze 
    130138      USE trc,               ONLY: ln_diatrc, med_diag, trc2d, trc3d  
     139      USE trcoxy_medusa,     ONLY: oxy_sato 
    131140 
    132141   !!* Substitution 
     
    140149      !! Loop avariables 
    141150      INTEGER :: ji, jj, jn 
    142  
    143       REAL(wp), DIMENSION(jpi,jpj) :: fregen,fregensi 
    144151 
    145152# if defined key_trc_diabio 
     
    158165      !! 
    159166      IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio' 
     167# endif 
     168 
     169      !! The section below, down to calculation of zo2min, was moved  
     170      !! from before the call to AIR_SEA in trcbio_medusa.F90 - marc 9/5/17  
     171      IF( lk_iomput ) THEN 
     172         DO jj = 2,jpjm1 
     173            DO ji = 2,jpim1 
     174               if (tmask(ji,jj,1) == 1) then 
     175                  !! sum tracers for inventory checks 
     176                  IF ( med_diag%INVTN%dgsave )   THEN 
     177                     ftot_n(ji,jj)  = ftot_n(ji,jj) +                        & 
     178                        (fse3t(ji,jj,jk) * (zphn(ji,jj) + zphd(ji,jj) +      & 
     179                                            zzmi(ji,jj) + zzme(ji,jj) +      & 
     180                                            zdet(ji,jj) + zdin(ji,jj))) 
     181                  ENDIF 
     182                  IF ( med_diag%INVTSI%dgsave )  THEN 
     183                     ftot_si(ji,jj) = ftot_si(ji,jj) +                       &  
     184                       (fse3t(ji,jj,jk) * (zpds(ji,jj) + zsil(ji,jj))) 
     185                  ENDIF 
     186                  IF ( med_diag%INVTFE%dgsave )  THEN 
     187                     ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       &  
     188                        (fse3t(ji,jj,jk) * (xrfn *                           & 
     189                                            (zphn(ji,jj) + zphd(ji,jj) +     & 
     190                                             zzmi(ji,jj) + zzme(ji,jj) +     & 
     191                                             zdet(ji,jj)) +                  & 
     192                                            zfer(ji,jj))) 
     193                  ENDIF 
     194               ENDIF 
     195            ENDDO 
     196         ENDDO 
     197 
     198# if defined key_roam 
     199         DO jj = 2,jpjm1 
     200            DO ji = 2,jpim1 
     201               if (tmask(ji,jj,1) == 1) then 
     202                  IF ( med_diag%INVTC%dgsave )  THEN 
     203                     ftot_c(ji,jj)  = ftot_c(ji,jj) +                        &  
     204                        (fse3t(ji,jj,jk) * ((xthetapn * zphn(ji,jj)) +       & 
     205                                            (xthetapd * zphd(ji,jj)) +       & 
     206                                            (xthetazmi * zzmi(ji,jj)) +      & 
     207                                            (xthetazme * zzme(ji,jj)) +      & 
     208                                            zdtc(ji,jj) + zdic(ji,jj))) 
     209                  ENDIF 
     210                  IF ( med_diag%INVTALK%dgsave ) THEN 
     211                     ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     & 
     212                                                       zalk(ji,jj)) 
     213                  ENDIF 
     214                  IF ( med_diag%INVTO2%dgsave )  THEN 
     215                     ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    & 
     216                                                        zoxy(ji,jj)) 
     217                  ENDIF 
     218               ENDIF 
     219            ENDDO 
     220         ENDDO 
     221 
     222         DO jj = 2,jpjm1 
     223            DO ji = 2,jpim1 
     224               if (tmask(ji,jj,1) == 1) then 
     225                  IF ( med_diag%INVTC%dgsave )  THEN 
     226                     !! 
     227                     !! AXY (10/11/16): CMIP6 diagnostics 
     228                     IF ( med_diag%INTDISSIC%dgsave ) THEN 
     229                        intdissic(ji,jj) = intdissic(ji,jj) +                & 
     230                                           (fse3t(ji,jj,jk) * zdic(ji,jj)) 
     231                     ENDIF 
     232                     IF ( med_diag%INTDISSIN%dgsave ) THEN 
     233                        intdissin(ji,jj) = intdissin(ji,jj) +                & 
     234                                           (fse3t(ji,jj,jk) * zdin(ji,jj)) 
     235                     ENDIF 
     236                     IF ( med_diag%INTDISSISI%dgsave ) THEN 
     237                        intdissisi(ji,jj) = intdissisi(ji,jj) +              & 
     238                                            (fse3t(ji,jj,jk) * zsil(ji,jj)) 
     239                     ENDIF 
     240                     IF ( med_diag%INTTALK%dgsave ) THEN 
     241                        inttalk(ji,jj) = inttalk(ji,jj) +                    & 
     242                                         (fse3t(ji,jj,jk) * zalk(ji,jj)) 
     243                     ENDIF 
     244                  ENDIF 
     245               ENDIF 
     246            ENDDO 
     247         ENDDO 
     248 
     249         DO jj = 2,jpjm1 
     250            DO ji = 2,jpim1 
     251               if (tmask(ji,jj,1) == 1) then 
     252                  IF ( med_diag%O2min%dgsave ) THEN 
     253                     if ( zoxy(ji,jj) < o2min(ji,jj) ) then 
     254                        o2min(ji,jj)  = zoxy(ji,jj) 
     255                        IF ( med_diag%ZO2min%dgsave ) THEN 
     256                           !! layer midpoint 
     257                           zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               & 
     258                                            fdep1(ji,jj)) / 2.0 
     259                        ENDIF 
     260                     endif 
     261                  ENDIF 
     262               ENDIF 
     263            ENDDO 
     264         ENDDO 
     265# endif 
     266      ENDIF 
     267 
     268# if defined key_roam 
     269      !! This section is moved from just below CALL to AIR_SEA in 
     270      !! trcbio_medusa.F90 - marc 9/5/17 
     271      DO jj = 2,jpjm1 
     272         DO ji = 2,jpim1 
     273            !! OPEN wet point IF..THEN loop 
     274            if (tmask(ji,jj,jk) == 1) then 
     275 
     276               !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 
     277               IF ( med_diag%O2SAT3%dgsave ) THEN 
     278! Remove f_o2sat3 - marc 9/5/17 
     279!                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 ) 
     280!                  o2sat3(ji, jj, jk) = f_o2sat3 
     281                  call oxy_sato( ztmp(ji,jj), zsal(ji,jj),                   & 
     282                                 o2sat3(ji,jj,jk) ) 
     283               ENDIF 
     284            ENDIF 
     285         ENDDO 
     286      ENDDO 
    160287# endif 
    161288 
  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90

    r7996 r8012  
    8585   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fslown, fslowc 
    8686   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fslownflux, fslowcflux 
     87   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregen,fregensi 
    8788   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregenfast,fregenfastsi 
    8889# if defined key_roam 
     
    309310               fdpn2(jpi,jpj),fdpd2(jpi,jpj),fdpds2(jpi,jpj),         & 
    310311               fdzmi2(jpi,jpj),fdzme2(jpi,jpj),                       & 
    311                fslown(jpi,jpj), fslowc(jpi,jpj),                      & 
     312               fslown(jpi,jpj),fslowc(jpi,jpj),                       & 
    312313               fslownflux(jpi,jpj),fslowcflux(jpi,jpj),               & 
     314               fregen(jpi,jpj),fregensi(jpi,jpj),                     & 
    313315               fregenfast(jpi,jpj),fregenfastsi(jpi,jpj),             & 
    314316# if defined key_roam 
  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90

    r7986 r8012  
    3737                                   f_fbenin_n, f_fbenin_si,                & 
    3838                                   f_omcal, fcaco3,                        & 
    39                                    fccd, fccd_dep, fdep1,                  & 
    40                                    fdpd, fdpds, fdzme, fgmepds,            & 
     39                                   fccd, fccd_dep, fdep1, fdd,             & 
     40                                   fdpd, fdpd2, fdpds, fdpds2,             & 
     41                                   fdpn, fdpn2,                            & 
     42                                   fdzme, fdzme2, fdzmi, fdzmi2,           & 
    4143                                   ffast2slowc, ffast2slown,               & 
    4244                                   ffastc, ffastca, ffastfe, ffastn,       & 
    4345                                   ffastsi,                                & 
     46                                   fgmed, fgmepd, fgmepds, fgmepn,         & 
     47                                   fgmezmi,                                & 
     48                                   fgmid, fgmipn,                          & 
     49                                   ficme, ficmi,                           & 
    4450                                   fifd_c, fifd_fe, fifd_n, fifd_si,       & 
     51                                   finme, finmi,                           & 
     52                                   fmeexcr, fmiexcr,                       & 
    4553                                   fofd_fe, fofd_n, fofd_si,               & 
    4654                                   fprotf,                                 & 
    47                                    fregenfast, fregenfastsi,               & 
     55                                   fregen, fregenfast, fregenfastsi,       & 
     56                                   fregensi,                               & 
    4857                                   freminc, freminca, freminfe,            & 
    4958                                   freminn, freminsi,                      & 
     59                                   fsdiss,                                 & 
    5060                                   fsedc, fsedca, fsedn, fsedfe, fsedsi,   & 
    51                                    fslowc, fslown,                         & 
     61                                   fslowc, fslowcflux,                     & 
     62                                   fslown, fslownflux,                     & 
    5263                                   ftempc, ftempca, ftempfe, ftempn,       & 
    5364                                   ftempsi,                                & 
     
    5566                                   fofd_c, fregenfastc,                    & 
    5667# endif 
    57                                    idf, idfval 
     68                                   idf, idfval,                            & 
     69                                   zdet, zdtc 
    5870      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,         & 
    5971                                   gphit, mbathy, tmask 
     
    6577                                   jexport, jfdfate, jinorgben, jocalccd,  & 
    6678                                   jorgben, jp_tem, jrratio,               & 
    67                                    ocal_ccd, xcaco3a, xcaco3b,             & 
     79                                   ocal_ccd, vsed,                         & 
     80                                   xbetac, xbetan,                         & 
     81                                   xcaco3a, xcaco3b,                       & 
    6882                                   xfastc, xfastca, xfastsi,               & 
    6983                                   xfdfrac1, xfdfrac2, xfdfrac3,           & 
    7084                                   xmassc, xmassca, xmasssi,               & 
    71                                    xprotca, xprotsi, xrfn, xridg_r0,       & 
     85                                   xphi, xprotca, xprotsi,                 & 
     86                                   xrfn, xridg_r0,                         & 
    7287                                   xsedc, xsedca, xsedfe,xsedn, xsedsi,    & 
    73                                    xthetapd, xthetazme,                    & 
     88                                   xthetapd, xthetapn,                     & 
     89                                   xthetazme, xthetazmi,                   & 
    7490                                   zn_sed_c, zn_sed_ca, zn_sed_fe,         & 
    7591                                   zn_sed_n, zn_sed_si 
     
    89105      !! temporary variables 
    90106      REAL(wp) :: fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8 
     107 
     108      !! The two sections below, slow detritus creation and Nutrient  
     109      !! regeneration are moved from just above the CALL to DETRITUS 
     110      !! in trcbio_medusa.F90. 
     111      !!--------------------------------------------------------- 
     112      !! Slow detritus creation 
     113      !!--------------------------------------------------------- 
     114      DO jj = 2,jpjm1 
     115         DO ji = 2,jpim1 
     116            IF (tmask(ji,jj,1) == 1) THEN 
     117               !! 
     118               !! this variable integrates the creation of slow sinking  
     119               !! detritus to allow the split between fast and slow  
     120               !! detritus to be diagnosed 
     121               fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) +                 & 
     122                                ((1.0 - xfdfrac1) * fdpd(ji,jj)) +           & 
     123                                ((1.0 - xfdfrac2) * fdzme(ji,jj)) +          & 
     124                                ((1.0 - xbetan) *                            & 
     125                                 (finmi(ji,jj) + finme(ji,jj))) 
     126               !! 
     127               !! this variable records the slow detrital sinking flux at  
     128               !! this particular depth; it is used in the output of this  
     129               !! flux at standard depths in the diagnostic outputs;  
     130               !! needs to be adjusted from per second to per day because  
     131               !! of parameter vsed 
     132               fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400. 
     133# if defined key_roam 
     134               !! 
     135               !! and the same for detrital carbon 
     136               fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) +                   & 
     137                                (xthetazmi * fdzmi(ji,jj)) +                 & 
     138                                (xthetapd * (1.0 - xfdfrac1) *               & 
     139                                 fdpd(ji,jj)) +                              & 
     140                                (xthetazme * (1.0 - xfdfrac2) *              & 
     141                                 fdzme(ji,jj)) +                             & 
     142                                ((1.0 - xbetac) * (ficmi(ji,jj) +            & 
     143                                                   ficme(ji,jj))) 
     144               !! 
     145               !! this variable records the slow detrital sinking flux  
     146               !! at this particular depth; it is used in the output of  
     147               !! this flux at standard depths in the diagnostic  
     148               !! outputs; needs to be adjusted from per second to per  
     149               !! day because of parameter vsed 
     150               fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400. 
     151# endif 
     152            ENDIF 
     153         ENDDO 
     154      ENDDO 
     155 
     156      !!--------------------------------------------------------- 
     157      !! Nutrient regeneration 
     158      !! this variable integrates total nitrogen regeneration down the 
     159      !! watercolumn; its value is stored and output as a 2D diagnostic; 
     160      !! the corresponding dissolution flux of silicon (from sources 
     161      !! other than fast detritus) is also integrated; note that, 
     162      !! confusingly, the linear loss terms from plankton compartments 
     163      !! are labelled as fdX2 when one might have expected fdX or fdX1 
     164      !!--------------------------------------------------------- 
     165      DO jj = 2,jpjm1 
     166         DO ji = 2,jpim1 
     167            IF (tmask(ji,jj,1) == 1) THEN 
     168               !! 
     169               !! nitrogen 
     170               fregen(ji,jj) =                                             & 
     171                                     ! messy feeding 
     172                               (((xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) + & 
     173                                 (xphi *                                   & 
     174                                  (fgmepn(ji,jj) + fgmepd(ji,jj) +         & 
     175                                   fgmezmi(ji,jj) + fgmed(ji,jj))) +       & 
     176                                     ! excretion + D remin. 
     177                                 fmiexcr(ji,jj) + fmeexcr(ji,jj) +         & 
     178                                 fdd(ji,jj) +                              & 
     179                                     ! linear mortality 
     180                                 fdpn2(ji,jj) + fdpd2(ji,jj) +             & 
     181                                 fdzmi2(ji,jj) + fdzme2(ji,jj)) *          & 
     182                                fse3t(ji,jj,jk)) 
     183               !! 
     184               !! silicon 
     185               fregensi(ji,jj) =                                           & 
     186                                     ! dissolution + non-lin. mortality 
     187                                 ((fsdiss(ji,jj) +                         & 
     188                                   ((1.0 - xfdfrac1) * fdpds(ji,jj)) +     & 
     189                                     ! egestion by zooplankton 
     190                                   ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +   & 
     191                                     ! linear mortality 
     192                                   fdpds2(ji,jj)) * fse3t(ji,jj,jk)) 
     193            ENDIF 
     194         ENDDO 
     195      ENDDO 
    91196 
    92197      !!------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90

    r7975 r8012  
    3535      USE bio_medusa_mod,    ONLY: fdpd, fdpd2, fdpds, fdpds2,             & 
    3636                                   fdpn, fdpn2, fdzme, fdzme2,             & 
    37                                    fdzmi, fdzmi2, fsin,                    & 
     37                                   fdzmi, fdzmi2, fsdiss, fsin,            & 
    3838                                   zphd, zphn, zpds, zzme, zzmi 
    3939      USE dom_oce,           ONLY: tmask 
     
    4343                                   xkphd, xkphn, xkzme, xkzmi,             & 
    4444                                   xmetapd, xmetapn, xmetazme, xmetazmi,   & 
    45                                    xmpd, xmpn, xmzme, xmzmi 
     45                                   xmpd, xmpn, xmzme, xmzmi, xsdiss 
    4646      USE zooplankton_mod,   ONLY: zooplankton 
    4747 
     
    178178      ENDDO 
    179179 
     180      !! diatom frustule dissolution. This section is moved from just 
     181      !! below CALL to iron_chem_scav in trcbio_medusa.F90 - marc 9/5/17 
     182      DO jj = 2,jpjm1 
     183         DO ji = 2,jpim1 
     184            IF (tmask(ji,jj,jk) == 1) THEN 
     185               fsdiss(ji,jj)  = xsdiss * zpds(ji,jj) 
     186            ENDIF 
     187         ENDDO 
     188      ENDDO 
     189 
    180190   END SUBROUTINE plankton 
    181191 
  • branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r7998 r8012  
    243243!      REAL(wp), DIMENSION(jpi,jpj) ::    fslown, fslowc 
    244244!      REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux 
    245       REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi 
     245!      REAL(wp), DIMENSION(jpi,jpj) ::    fregen,fregensi 
    246246!      REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi 
    247247# if defined key_roam 
     
    397397      ENDIF 
    398398 
    399       !!---------------------------------------------------------------------- 
     399      !!------------------------------------------------------------------ 
    400400      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency 
    401401      !! terms of all biological equations to 0. 
    402       !!---------------------------------------------------------------------- 
     402      !!------------------------------------------------------------------ 
    403403      !! 
    404404      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit 
     
    410410      b0 = 1. 
    411411# endif 
    412       !!---------------------------------------------------------------------- 
     412      !!------------------------------------------------------------------ 
    413413      !! fast detritus ballast scheme (0 = no; 1 = yes) 
    414414      !! alternative to ballast scheme is same scheme but with no ballast 
    415415      !! protection (not dissimilar to Martin et al., 1987) 
    416       !!---------------------------------------------------------------------- 
     416      !!------------------------------------------------------------------ 
    417417      !! 
    418418      iball = 1 
    419419 
    420       !!---------------------------------------------------------------------- 
     420      !!------------------------------------------------------------------ 
    421421      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output 
    422422      !! these should *only* be used in 1D since they give comprehensive 
    423423      !! output for ecological functions in the model; primarily used in 
    424424      !! debugging 
    425       !!---------------------------------------------------------------------- 
     425      !!------------------------------------------------------------------ 
    426426      !! 
    427427      idf    = 0 
     
    434434      endif 
    435435 
    436       !!---------------------------------------------------------------------- 
     436      !!------------------------------------------------------------------ 
    437437      !! Initialise arrays to zero and set up arrays for diagnostics 
    438       !!---------------------------------------------------------------------- 
    439 ! tmp - marc 
    440       write(numout,*) 'bbb13. before call to bio_medusa_init,kt=',kt 
    441       flush(numout) 
    442 ! 
     438      !!------------------------------------------------------------------ 
    443439      CALL bio_medusa_init( kt ) 
    444 ! tmp - marc 
    445       write(numout,*) 'bbb14. after call to bio_medusa_init,kt=',kt 
    446       flush(numout) 
    447 ! 
    448        !! 
     440 
    449441# if defined key_axy_nancheck 
    450442       DO jn = 1,jptra 
     
    452444          !! fq1 = MAXVAL(trn(:,:,:,jn)) 
    453445          fq2 = SUM(trn(:,:,:,jn)) 
    454           !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', & 
    455           !! &        kt, jn, fq0, fq1, fq2 
    456           !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR 
    457           !!                 and has been replaced here with a specialist routine 
     446          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     & 
     447          !!                kt, jn, fq0, fq1, fq2 
     448          !! AXY (30/01/14): much to our surprise, the next line doesn't  
     449          !!                 work on HECTOR and has been replaced here with  
     450          !!                 a specialist routine 
    458451          !! if (fq2 /= fq2 ) then 
    459452          if ( ieee_is_nan( fq2 ) ) then 
    460453             !! there's a NaN here 
    461              if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:' 
     454             if (lwp) write(numout,*) 'NAN detected in field', jn,           & 
     455                                      'at time', kt, 'at position:' 
    462456             DO jk = 1,jpk 
    463457                DO jj = 1,jpj 
     
    466460                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then 
    467461                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then 
    468                          if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', & 
    469                          &        tmask(ji,jj,jk), ji, jj, jk, jn 
     462                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           & 
     463                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn 
    470464                      endif 
    471465                   enddo 
     
    479473 
    480474# if defined key_debug_medusa 
    481       IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 
     475      IF (lwp) write (numout,*)                                              & 
     476                     'trc_bio_medusa: variables initialised and checked' 
    482477      CALL flush(numout) 
    483478# endif  
    484479 
    485480# if defined key_roam 
    486       !!---------------------------------------------------------------------- 
     481      !!------------------------------------------------------------------ 
    487482      !! calculate atmospheric pCO2 
    488       !!---------------------------------------------------------------------- 
     483      !!------------------------------------------------------------------ 
    489484      !! 
    490485      !! what's atmospheric pCO2 doing? (data start in 1859) 
     
    504499         !! AXY (14/06/12): tweaked to make more sense (and be correct) 
    505500#  if defined key_bs_axy_yrlen 
    506          fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  !! bugfix: for 360d year with HadGEM2-ES forcing 
     501         !! bugfix: for 360d year with HadGEM2-ES forcing 
     502         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0   
    507503#  else 
    508          fq3 = (real(nday_year) - 1.0 + fq2) / 365.0  !! original use of 365 days (not accounting for leap year or 360d year) 
     504         !! original use of 365 days (not accounting for leap year or  
     505         !! 360d year) 
     506         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 
    509507#  endif 
    510508         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
     
    512510      endif 
    513511#  if defined key_axy_pi_co2 
    514       f_xco2a(:,:) = 284.725          !! OCMIP pre-industrial pCO2 
     512      !! OCMIP pre-industrial pCO2 
     513      f_xco2a(:,:) = 284.725 
    515514#  endif 
    516515      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
     
    540539      !!============================= 
    541540      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call : 
    542       !!          we don't want to call on the first time-step of all run submission,  
    543       !!          but only on the very first time-step, and then every month 
    544       !!          So we call on nittrc000 if not restarted run,  
    545       !!          else if one month after last call. 
    546       !!          assume one month is 30d --> 3600*24*30 : 2592000s 
    547       !!          try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)    
     541      !!          we don't want to call on the first time-step of all run  
     542      !!          submission, but only on the very first time-step, and  
     543      !!          then every month. So we call on nittrc000 if not  
     544      !!          restarted run, else if one month after last call. 
     545      !!          Assume one month is 30d --> 3600*24*30 : 2592000s 
     546      !!          try to call carb-chem at 1st month's tm-stp :  
     547      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    548548      !!          ++ need to pass carb-chem output var through restarts 
    549       If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN 
    550          !!---------------------------------------------------------------------- 
     549      If ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR.                      & 
     550           ( mod(kt*rdt,2592000.) == rdt ) ) THEN 
     551         !!--------------------------------------------------------------- 
    551552         !! Calculate the carbonate chemistry for the whole ocean on the first 
    552553         !! simulation timestep and every month subsequently; the resulting 3D 
    553554         !! field of omega calcite is used to determine the depth of the CCD 
    554          !!---------------------------------------------------------------------- 
     555         !!--------------------------------------------------------------- 
    555556         CALL carb_chem( kt ) 
    556557 
     
    563564# endif  
    564565 
    565       !!---------------------------------------------------------------------- 
     566      !!------------------------------------------------------------------ 
    566567      !! MEDUSA has unified equation through the water column 
    567568      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)  
    568569      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1           
    569       !!---------------------------------------------------------------------- 
     570      !!------------------------------------------------------------------ 
    570571      !! 
    571572      !! NOTE: the ordering of the loops below differs from that of some other 
     
    583584         !! OPEN horizontal loops 
    584585         DO jj = 2,jpjm1 
    585          DO ji = 2,jpim1 
    586             !! OPEN wet point IF..THEN loop 
    587             if (tmask(ji,jj,jk) == 1) then                
    588                !!=========================================================== 
    589                !! SETUP LOCAL GRID CELL 
    590                !!=========================================================== 
    591                !! 
    592                !!----------------------------------------------------------- 
    593                !! Some notes on grid vertical structure 
    594                !! - fsdepw(ji,jj,jk) is the depth of the upper surface of  
    595                !!   level jk 
    596                !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of  
    597                !!   level jk 
    598                !! - fse3t(ji,jj,jk)  is the thickness of level jk 
    599                !!----------------------------------------------------------- 
    600                !! 
    601                !! AXY (01/03/10): set up level depth (bottom of level) 
    602                fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 
    603                !! AXY (28/11/16): local seafloor depth 
    604                !!                 previously mbathy(ji,jj) - 1, now  
    605                !!                 mbathy(ji,jj) 
    606 ! I should be able to remove this - marc 5/5/17 
    607 !               mbathy(ji,jj) = mbathy(ji,jj) 
    608                !! 
    609                !! set up model tracers 
    610                !! negative values of state variables are not allowed to 
    611                !! contribute to the calculated fluxes 
    612                !! non-diatom chlorophyll 
    613                zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 
    614                !! diatom chlorophyll 
    615                zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 
    616                !! non-diatoms 
    617                zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 
    618                !! diatoms 
    619                zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 
    620                !! diatom silicon 
    621                zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 
    622                !! AXY (28/01/10): probably need to take account of  
    623                !! chl/biomass connection 
    624                if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 
    625                if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. 
    626                if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0. 
    627                if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0. 
    628           !! AXY (23/01/14): duh - why did I forget diatom silicon? 
    629           if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 
    630           if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 
    631                !! microzooplankton 
    632                zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 
    633                !! mesozooplankton 
    634                zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 
    635                !! detrital nitrogen 
    636                zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 
    637                !! dissolved inorganic nitrogen 
    638                zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 
    639                !! dissolved silicic acid 
    640                zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 
    641                !! dissolved "iron" 
    642                zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 
     586            DO ji = 2,jpim1 
     587               !! OPEN wet point IF..THEN loop 
     588               if (tmask(ji,jj,jk) == 1) then                
     589                  !!========================================================= 
     590                  !! SETUP LOCAL GRID CELL 
     591                  !!========================================================= 
     592                  !! 
     593                  !!--------------------------------------------------------- 
     594                  !! Some notes on grid vertical structure 
     595                  !! - fsdepw(ji,jj,jk) is the depth of the upper surface of  
     596                  !!   level jk 
     597                  !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of  
     598                  !!   level jk 
     599                  !! - fse3t(ji,jj,jk)  is the thickness of level jk 
     600                  !!--------------------------------------------------------- 
     601                  !! 
     602                  !! AXY (01/03/10): set up level depth (bottom of level) 
     603                  fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 
     604                  !! 
     605                  !! set up model tracers 
     606                  !! negative values of state variables are not allowed to 
     607                  !! contribute to the calculated fluxes 
     608                  !! non-diatom chlorophyll 
     609                  zchn(ji,jj) = max(0.,trn(ji,jj,jk,jpchn)) 
     610                  !! diatom chlorophyll 
     611                  zchd(ji,jj) = max(0.,trn(ji,jj,jk,jpchd)) 
     612                  !! non-diatoms 
     613                  zphn(ji,jj) = max(0.,trn(ji,jj,jk,jpphn)) 
     614                  !! diatoms 
     615                  zphd(ji,jj) = max(0.,trn(ji,jj,jk,jpphd)) 
     616                  !! diatom silicon 
     617                  zpds(ji,jj) = max(0.,trn(ji,jj,jk,jppds)) 
     618                  !! AXY (28/01/10): probably need to take account of  
     619                  !! chl/biomass connection 
     620                  if (zchn(ji,jj).eq.0.) zphn(ji,jj) = 0. 
     621                  if (zchd(ji,jj).eq.0.) zphd(ji,jj) = 0. 
     622                  if (zphn(ji,jj).eq.0.) zchn(ji,jj) = 0. 
     623                  if (zphd(ji,jj).eq.0.) zchd(ji,jj) = 0. 
     624             !! AXY (23/01/14): duh - why did I forget diatom silicon? 
     625             if (zpds(ji,jj).eq.0.) zphd(ji,jj) = 0. 
     626             if (zphd(ji,jj).eq.0.) zpds(ji,jj) = 0. 
     627               ENDIF 
     628            ENDDO 
     629         ENDDO 
     630 
     631         DO jj = 2,jpjm1 
     632            DO ji = 2,jpim1 
     633               if (tmask(ji,jj,1) == 1) then 
     634                  !! microzooplankton 
     635                  zzmi(ji,jj) = max(0.,trn(ji,jj,jk,jpzmi)) 
     636                  !! mesozooplankton 
     637                  zzme(ji,jj) = max(0.,trn(ji,jj,jk,jpzme)) 
     638                  !! detrital nitrogen 
     639                  zdet(ji,jj) = max(0.,trn(ji,jj,jk,jpdet)) 
     640                  !! dissolved inorganic nitrogen 
     641                  zdin(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) 
     642                  !! dissolved silicic acid 
     643                  zsil(ji,jj) = max(0.,trn(ji,jj,jk,jpsil)) 
     644                  !! dissolved "iron" 
     645                  zfer(ji,jj) = max(0.,trn(ji,jj,jk,jpfer)) 
     646               ENDIF 
     647            ENDDO 
     648         ENDDO 
     649 
    643650# if defined key_roam 
    644                !! detrital carbon 
    645                zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 
    646                !! dissolved inorganic carbon 
    647                zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 
    648                !! alkalinity 
    649                zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 
    650                !! oxygen 
    651                zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 
     651         DO jj = 2,jpjm1 
     652            DO ji = 2,jpim1 
     653               if (tmask(ji,jj,1) == 1) then 
     654                  !! detrital carbon 
     655                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 
     656                  !! dissolved inorganic carbon 
     657                  zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 
     658                  !! alkalinity 
     659                  zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 
     660                  !! oxygen 
     661                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 
    652662#  if defined key_axy_carbchem && defined key_mocsy 
    653                !! phosphate via DIN and Redfield 
    654                zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
     663                  !! phosphate via DIN and Redfield 
     664                  zpho(ji,jj) = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 
    655665#  endif 
    656                !! 
    657                !! also need physical parameters for gas exchange calculations 
    658                ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 
    659                zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 
    660                !! 
    661           !! AXY (28/02/14): check input fields 
    662                if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 
    663                   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ',  & 
    664                      tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',     & 
    665                      ji, ',', jj, ',', jk, ') at time', kt 
    666         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ',& 
    667                      tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    668                   !! temperatur 
    669                   ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 
    670                endif 
    671                if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 
    672                   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 
    673                      tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    & 
    674                      ji, ',', jj, ',', jk, ') at time', kt 
    675                endif 
     666                  !! 
     667                  !! also need physical parameters for gas exchange  
     668                  !! calculations 
     669                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 
     670                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 
     671                  !! 
     672             !! AXY (28/02/14): check input fields 
     673                  if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 
     674                     IF(lwp) WRITE(numout,*)                                 & 
     675                        ' trc_bio_medusa: T WARNING 2D, ',                   & 
     676                        tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          & 
     677                        ' at (', ji, ',', jj, ',', jk, ') at time', kt 
     678           IF(lwp) WRITE(numout,*)                                 & 
     679                        ' trc_bio_medusa: T SWITCHING 2D, ',                 & 
     680                        tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
     681                     !! temperatur 
     682                     ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 
     683                  endif 
     684                  if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 
     685                     IF(lwp) WRITE(numout,*)                                 & 
     686                        ' trc_bio_medusa: S WARNING 2D, ',                   & 
     687                        tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          & 
     688                        ' at (', ji, ',', jj, ',', jk, ') at time', kt 
     689                  endif 
     690               ENDIF 
     691            ENDDO 
     692         ENDDO 
    676693# else 
    677                !! implicit detrital carbon 
    678                zdtc(ji,jj) = zdet(ji,jj) * xthetad 
     694         DO jj = 2,jpjm1 
     695            DO ji = 2,jpim1 
     696               if (tmask(ji,jj,1) == 1) then 
     697                  !! implicit detrital carbon 
     698                  zdtc(ji,jj) = zdet(ji,jj) * xthetad 
     699               ENDIF 
     700            ENDDO 
     701         ENDDO 
    679702# endif 
    680703# if defined key_debug_medusa 
    681                if (idf.eq.1) then 
    682                   !! AXY (15/01/10) 
    683                   if (trn(ji,jj,jk,jpdin).lt.0.) then 
    684                      IF (lwp) write (numout,*) '------------------------------' 
    685                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',        & 
    686                         trn(ji,jj,jk,jpdin) 
    687                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',        & 
    688                         ji, jj, jk, kt 
     704         DO jj = 2,jpjm1 
     705            DO ji = 2,jpim1 
     706               if (tmask(ji,jj,1) == 1) then 
     707                  if (idf.eq.1) then 
     708                     !! AXY (15/01/10) 
     709                     if (trn(ji,jj,jk,jpdin).lt.0.) then 
     710                        IF (lwp) write (numout,*)                            & 
     711                           '------------------------------' 
     712                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =',    & 
     713                           trn(ji,jj,jk,jpdin) 
     714                        IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @',    & 
     715                           ji, jj, jk, kt 
     716                     endif 
     717                     if (trn(ji,jj,jk,jpsil).lt.0.) then 
     718                        IF (lwp) write (numout,*)                            & 
     719                           '------------------------------' 
     720                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',    & 
     721                           trn(ji,jj,jk,jpsil) 
     722                        IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',    & 
     723                           ji, jj, jk, kt 
     724                     endif 
     725#  if defined key_roam 
     726                     if (trn(ji,jj,jk,jpdic).lt.0.) then 
     727                        IF (lwp) write (numout,*)                            & 
     728                           '------------------------------' 
     729                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',    & 
     730                           trn(ji,jj,jk,jpdic) 
     731                        IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',    & 
     732                           ji, jj, jk, kt 
     733                     endif 
     734                     if (trn(ji,jj,jk,jpalk).lt.0.) then 
     735                        IF (lwp) write (numout,*)                            & 
     736                           '------------------------------' 
     737                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',    & 
     738                           trn(ji,jj,jk,jpalk) 
     739                        IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',    & 
     740                           ji, jj, jk, kt 
     741                     endif 
     742                     if (trn(ji,jj,jk,jpoxy).lt.0.) then 
     743                        IF (lwp) write (numout,*)                            & 
     744                           '------------------------------' 
     745                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',    & 
     746                           trn(ji,jj,jk,jpoxy) 
     747                        IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',    & 
     748                           ji, jj, jk, kt 
     749                     endif 
     750#  endif 
    689751                  endif 
    690                   if (trn(ji,jj,jk,jpsil).lt.0.) then 
    691                      IF (lwp) write (numout,*) '------------------------------' 
    692                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =',        & 
    693                         trn(ji,jj,jk,jpsil) 
    694                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @',        & 
    695                         ji, jj, jk, kt 
    696                   endif 
     752               ENDIF 
     753            ENDDO 
     754         ENDDO 
     755# endif 
     756# if defined key_debug_medusa 
     757! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17 
     758!         if (idf.eq.1.AND.idfval.eq.1) then 
     759!            DO jj = 2,jpjm1 
     760!               DO ji = 2,jpim1 
     761!                  if (tmask(ji,jj,1) == 1) then 
     762!                     !! report state variable values 
     763!                     IF (lwp) write (numout,*)                               & 
     764!                        '------------------------------' 
     765!                     IF (lwp) write (numout,*) 'fthk(',jk,') = ',            & 
     766!                        fse3t(ji,jj,jk) 
     767!                     IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj) 
     768!                     IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj) 
     769!                     IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj) 
     770!                     IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj) 
     771!                     IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj) 
     772!                     IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj) 
     773!                     IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj) 
     774!                     IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj) 
     775!                     IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj) 
    697776#  if defined key_roam 
    698                   if (trn(ji,jj,jk,jpdic).lt.0.) then 
    699                      IF (lwp) write (numout,*) '------------------------------' 
    700                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =',        & 
    701                         trn(ji,jj,jk,jpdic) 
    702                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @',        & 
    703                         ji, jj, jk, kt 
    704                   endif 
    705                   if (trn(ji,jj,jk,jpalk).lt.0.) then 
    706                      IF (lwp) write (numout,*) '------------------------------' 
    707                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =',        & 
    708                         trn(ji,jj,jk,jpalk) 
    709                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @',        & 
    710                         ji, jj, jk, kt 
    711                   endif 
    712                   if (trn(ji,jj,jk,jpoxy).lt.0.) then 
    713                      IF (lwp) write (numout,*) '------------------------------' 
    714                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =',        & 
    715                         trn(ji,jj,jk,jpoxy) 
    716                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @',        & 
    717                         ji, jj, jk, kt 
    718                   endif 
     777!                     IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj) 
     778!                     IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj) 
     779!                     IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj) 
     780!                     IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj) 
    719781#  endif 
    720                endif 
    721 # endif 
     782!                  ENDIF 
     783!               ENDDO 
     784!            ENDDO 
     785!         endif 
     786# endif 
     787 
    722788# if defined key_debug_medusa 
    723                !! report state variable values 
    724                if (idf.eq.1.AND.idfval.eq.1) then 
    725                   IF (lwp) write (numout,*) '------------------------------' 
    726                   IF (lwp) write (numout,*) 'fthk(',jk,') = ', fse3t(ji,jj,jk) 
    727                   IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn(ji,jj) 
    728                   IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd(ji,jj) 
    729                   IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds(ji,jj) 
    730                   IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi(ji,jj) 
    731                   IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme(ji,jj) 
    732                   IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet(ji,jj) 
    733                   IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin(ji,jj) 
    734                   IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil(ji,jj) 
    735                   IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer(ji,jj) 
    736 #  if defined key_roam 
    737                   IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc(ji,jj) 
    738                   IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic(ji,jj) 
    739                   IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk(ji,jj) 
    740                   IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy(ji,jj)                   
    741 #  endif 
    742                endif 
    743 # endif 
    744  
    745 # if defined key_debug_medusa 
    746                if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 
    747                   IF (lwp) write (numout,*) '------------------------------' 
    748                   IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj) 
    749                endif 
    750 # endif 
    751  
    752                !! sum tracers for inventory checks 
    753                IF( lk_iomput ) THEN 
    754                   IF ( med_diag%INVTN%dgsave )   THEN 
    755                      ftot_n(ji,jj)  = ftot_n(ji,jj) +                        & 
    756                         (fse3t(ji,jj,jk) * ( zphn(ji,jj) + zphd(ji,jj) +     & 
    757                                              zzmi(ji,jj) + zzme(ji,jj) +     & 
    758                                              zdet(ji,jj) + zdin(ji,jj) ) ) 
    759                   ENDIF 
    760                   IF ( med_diag%INVTSI%dgsave )  THEN 
    761                      ftot_si(ji,jj) = ftot_si(ji,jj) +                       &  
    762                         (fse3t(ji,jj,jk) * ( zpds(ji,jj) + zsil(ji,jj) ) ) 
    763                   ENDIF 
    764                   IF ( med_diag%INVTFE%dgsave )  THEN 
    765                      ftot_fe(ji,jj) = ftot_fe(ji,jj) +                       &  
    766                         (fse3t(ji,jj,jk) * ( xrfn *                          & 
    767                                             ( zphn(ji,jj) + zphd(ji,jj) +    & 
    768                                               zzmi(ji,jj) + zzme(ji,jj) +    & 
    769                                               zdet(ji,jj) ) + zfer(ji,jj) ) ) 
    770                   ENDIF 
    771 # if defined key_roam 
    772                   IF ( med_diag%INVTC%dgsave )  THEN 
    773                      ftot_c(ji,jj)  = ftot_c(ji,jj) + &  
    774                         (fse3t(ji,jj,jk) * ( (xthetapn * zphn(ji,jj)) +      & 
    775                                              (xthetapd * zphd(ji,jj)) +      & 
    776                                              (xthetazmi * zzmi(ji,jj)) +     & 
    777                                              (xthetazme * zzme(ji,jj)) +     & 
    778                                              zdtc(ji,jj) + zdic(ji,jj) ) ) 
    779                   ENDIF 
    780                   IF ( med_diag%INVTALK%dgsave ) THEN 
    781                      ftot_a(ji,jj)  = ftot_a(ji,jj) + (fse3t(ji,jj,jk) *     & 
    782                                                        zalk(ji,jj)) 
    783                   ENDIF 
    784                   IF ( med_diag%INVTO2%dgsave )  THEN 
    785                      ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fse3t(ji,jj,jk) *    & 
    786                                                         zoxy(ji,jj)) 
    787                   ENDIF 
    788                   !! 
    789                   !! AXY (10/11/16): CMIP6 diagnostics 
    790                   IF ( med_diag%INTDISSIC%dgsave ) THEN 
    791                      intdissic(ji,jj) = intdissic(ji,jj) +                   & 
    792                                         (fse3t(ji,jj,jk) * zdic(ji,jj)) 
    793                   ENDIF 
    794                   IF ( med_diag%INTDISSIN%dgsave ) THEN 
    795                      intdissin(ji,jj) = intdissin(ji,jj) +                   & 
    796                                         (fse3t(ji,jj,jk) * zdin(ji,jj)) 
    797                   ENDIF 
    798                   IF ( med_diag%INTDISSISI%dgsave ) THEN 
    799                      intdissisi(ji,jj) = intdissisi(ji,jj) +                 & 
    800                                          (fse3t(ji,jj,jk) * zsil(ji,jj)) 
    801                   ENDIF 
    802                   IF ( med_diag%INTTALK%dgsave ) THEN 
    803                      inttalk(ji,jj) = inttalk(ji,jj) +                       & 
    804                                       (fse3t(ji,jj,jk) * zalk(ji,jj)) 
    805                   ENDIF 
    806                   IF ( med_diag%O2min%dgsave ) THEN 
    807                      if ( zoxy(ji,jj) < o2min(ji,jj) ) then 
    808                         o2min(ji,jj)  = zoxy(ji,jj) 
    809                         IF ( med_diag%ZO2min%dgsave ) THEN 
    810                            !! layer midpoint 
    811                            zo2min(ji,jj) = (fsdepw(ji,jj,jk) +               & 
    812                                             fdep1(ji,jj)) / 2.0 
    813                         ENDIF 
    814                      endif 
    815                   ENDIF 
    816 # endif 
    817                ENDIF 
    818  
    819                CALL flush(numout) 
    820  
    821             ENDIF 
    822          ENDDO 
    823          ENDDO 
    824  
    825          !!---------------------------------------------------------------- 
     789! I'M NOT SURE THIS USEFUL NOW THAT I'VE SPLIT THE DO LOOP - marc 8/5/17 
     790!         if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 
     791!            DO jj = 2,jpjm1 
     792!               DO ji = 2,jpim1 
     793!                  if (tmask(ji,jj,1) == 1) then 
     794!                     IF (lwp) write (numout,*)                               & 
     795!                       '------------------------------' 
     796!                     IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj) 
     797!                  ENDIF 
     798!               ENDDO 
     799!            ENDDO 
     800!         endif 
     801# endif 
     802 
     803         !!--------------------------------------------------------------- 
    826804         !! Calculate air-sea gas exchange and river inputs 
    827          !!---------------------------------------------------------------- 
     805         !!--------------------------------------------------------------- 
    828806         IF ( jk == 1 ) THEN 
    829             call air_sea( kt ) 
    830          END IF 
    831  
    832 # if defined key_roam 
    833  
    834 ! Moved from above - marc 21/4/17 
    835 ! I think this should be moved into diagnostics at bottom - it 
    836 ! doesn't like it is used anyway else - marc 21/4/17 
    837          DO jj = 2,jpjm1 
    838          DO ji = 2,jpim1 
    839             !! OPEN wet point IF..THEN loop 
    840             if (tmask(ji,jj,jk) == 1) then 
    841  
    842                !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 
    843                IF ( med_diag%O2SAT3%dgsave ) THEN 
    844                   call oxy_sato( ztmp(ji,jj), zsal(ji,jj), f_o2sat3 ) 
    845                   o2sat3(ji, jj, jk) = f_o2sat3 
    846                ENDIF 
    847             ENDIF 
    848          ENDDO 
    849          ENDDO 
    850 # endif 
    851  
    852          !!------------------------------------------------------------------ 
     807            CALL air_sea( kt ) 
     808         ENDIF 
     809 
     810         !!--------------------------------------------------------------- 
    853811         !! Phytoplankton growth, zooplankton grazing and miscellaneous 
    854812         !! plankton losses.  
    855          !!------------------------------------------------------------------ 
     813         !!--------------------------------------------------------------- 
    856814         CALL plankton( jk ) 
    857815 
    858          !!---------------------------------------------------------------- 
     816         !!--------------------------------------------------------------- 
    859817         !! Iron chemistry and scavenging 
    860          !!---------------------------------------------------------------- 
     818         !!--------------------------------------------------------------- 
    861819         CALL iron_chem_scav( jk ) 
    862820 
    863 ! Miscellaneous processes - marc 
    864  
    865          DO jj = 2,jpjm1 
    866          DO ji = 2,jpim1 
    867             !! OPEN wet point IF..THEN loop 
    868             if (tmask(ji,jj,jk) == 1) then 
    869  
    870                !!---------------------------------------------------------------------- 
    871                !! Miscellaneous 
    872                !!---------------------------------------------------------------------- 
    873                !! 
    874                !! diatom frustule dissolution 
    875                fsdiss(ji,jj)  = xsdiss * zpds(ji,jj) 
    876  
    877 # if defined key_debug_medusa 
    878                !! report miscellaneous calculations 
    879                if (idf.eq.1.AND.idfval.eq.1) then 
    880                   IF (lwp) write (numout,*) '------------------------------' 
    881                   IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss(ji,jj) 
    882                endif 
    883 # endif 
    884  
    885                !!---------------------------------------------------------------------- 
    886                !! Slow detritus creation 
    887                !!---------------------------------------------------------------------- 
    888                !! this variable integrates the creation of slow sinking detritus 
    889                !! to allow the split between fast and slow detritus to be  
    890                !! diagnosed 
    891                fslown(ji,jj)  = fdpn(ji,jj) + fdzmi(ji,jj) + ((1.0 - xfdfrac1) * fdpd(ji,jj)) + & 
    892                ((1.0 - xfdfrac2) * fdzme(ji,jj)) + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj))) 
    893                !! 
    894                !! this variable records the slow detrital sinking flux at this 
    895                !! particular depth; it is used in the output of this flux at 
    896                !! standard depths in the diagnostic outputs; needs to be 
    897                !! adjusted from per second to per day because of parameter vsed 
    898                fslownflux(ji,jj) = zdet(ji,jj) * vsed * 86400. 
    899 # if defined key_roam 
    900                !! 
    901                !! and the same for detrital carbon 
    902                fslowc(ji,jj)  = (xthetapn * fdpn(ji,jj)) + (xthetazmi * fdzmi(ji,jj)) + & 
    903                (xthetapd * (1.0 - xfdfrac1) * fdpd(ji,jj)) + & 
    904                (xthetazme * (1.0 - xfdfrac2) * fdzme(ji,jj)) + & 
    905                ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj))) 
    906                !! 
    907                !! this variable records the slow detrital sinking flux at this 
    908                !! particular depth; it is used in the output of this flux at 
    909                !! standard depths in the diagnostic outputs; needs to be 
    910                !! adjusted from per second to per day because of parameter vsed 
    911                fslowcflux(ji,jj) = zdtc(ji,jj) * vsed * 86400. 
    912 # endif 
    913  
    914                !!---------------------------------------------------------------------- 
    915                !! Nutrient regeneration 
    916                !! this variable integrates total nitrogen regeneration down the 
    917                !! watercolumn; its value is stored and output as a 2D diagnostic; 
    918                !! the corresponding dissolution flux of silicon (from sources 
    919                !! other than fast detritus) is also integrated; note that, 
    920                !! confusingly, the linear loss terms from plankton compartments 
    921                !! are labelled as fdX2 when one might have expected fdX or fdX1 
    922                !!---------------------------------------------------------------------- 
    923                !! 
    924                !! nitrogen 
    925                fregen(ji,jj)   = (( (xphi * (fgmipn(ji,jj) + fgmid(ji,jj))) +                &  ! messy feeding 
    926                (xphi * (fgmepn(ji,jj) + fgmepd(ji,jj) + fgmezmi(ji,jj) + fgmed(ji,jj))) +           &  ! messy feeding 
    927                fmiexcr(ji,jj) + fmeexcr(ji,jj) + fdd(ji,jj) +                                &  ! excretion + D remin. 
    928                fdpn2(ji,jj) + fdpd2(ji,jj) + fdzmi2(ji,jj) + fdzme2(ji,jj)) * fse3t(ji,jj,jk))                    ! linear mortality 
    929                !! 
    930                !! silicon 
    931                fregensi(ji,jj) = (( fsdiss(ji,jj) + ((1.0 - xfdfrac1) * fdpds(ji,jj)) +      &  ! dissolution + non-lin. mortality 
    932                ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +                           &  ! egestion by zooplankton 
    933                fdpds2(ji,jj)) * fse3t(ji,jj,jk))                                             ! linear mortality 
    934 # if defined key_roam 
    935                !! 
    936                !! carbon 
    937 ! Doesn't look this is used - marc 10/4/17 
    938 !               fregenc(ji,jj)  = (( (xphi * ((xthetapn * fgmipn(ji,jj)) + fgmidc(ji,jj))) +  &  ! messy feeding 
    939 !               (xphi * ((xthetapn * fgmepn(ji,jj)) + (xthetapd * fgmepd(ji,jj)) +     &  ! messy feeding 
    940 !               (xthetazmi * fgmezmi(ji,jj)) + fgmedc(ji,jj))) +                       &  ! messy feeding 
    941 !               fmiresp(ji,jj) + fmeresp(ji,jj) + fddc(ji,jj) +                               &  ! respiration + D remin. 
    942 !               (xthetapn * fdpn2(ji,jj)) + (xthetapd * fdpd2(ji,jj)) +                &  ! linear mortality 
    943 !               (xthetazmi * fdzmi2(ji,jj)) + (xthetazme * fdzme2(ji,jj))) * fse3t(ji,jj,jk))        ! linear mortality 
    944 # endif 
    945  
    946             ENDIF 
    947          ENDDO 
    948          ENDDO 
    949  
    950          !!------------------------------------------------------------------ 
    951          !! Detritus process 
    952          !!------------------------------------------------------------------ 
     821         !!--------------------------------------------------------------- 
     822         !! Detritus processes 
     823         !!--------------------------------------------------------------- 
    953824         CALL detritus( jk, iball ) 
    954825 
    955          !!------------------------------------------------------------------ 
     826         !!--------------------------------------------------------------- 
    956827         !! Updating tracers 
    957          !!------------------------------------------------------------------ 
     828         !!--------------------------------------------------------------- 
    958829         CALL bio_medusa_update( kt, jk ) 
    959830 
    960          !!------------------------------------------------------------------ 
     831         !!--------------------------------------------------------------- 
    961832         !! Diagnostics 
    962          !!------------------------------------------------------------------ 
     833         !!--------------------------------------------------------------- 
    963834         CALL bio_medusa_diag( kt, jk ) 
    964835 
     836         !!------------------------------------------------------- 
     837         !! 2d specific k level diags 
     838         !!------------------------------------------------------- 
    965839         IF( lk_iomput  .AND.  .NOT.  ln_diatrc  ) THEN 
    966  
    967             !!------------------------------------------------------- 
    968             !! 2d specific k level diags 
    969             !!------------------------------------------------------- 
    970840            CALL bio_medusa_diag_slice( jk ) 
    971  
    972841         ENDIF 
    973842 
     
    975844      ENDDO 
    976845 
    977       !!---------------------------------------------------------------------- 
     846      !!------------------------------------------------------------------ 
    978847      !! Final calculations for diagnostics 
    979       !!---------------------------------------------------------------------- 
     848      !!------------------------------------------------------------------ 
    980849      CALL bio_medusa_fin( kt ) 
    981850 
Note: See TracChangeset for help on using the changeset viewer.