Ignore:
Timestamp:
2017-09-14T17:42:14+02:00 (3 years ago)
Author:
marc
Message:

gmed:#346 Merging branches/NERC/dev_r5518_GO6_MEDUSA_conserv

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r8442 r8521  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Amend bethic reservoir updating 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_medusa 
     
    6263      REAL(wp) :: fq0,fq1,fq2,fq3 
    6364 
     65# if defined key_roam                      
     66      !!---------------------------------------------------------------------- 
     67      !! AXY (09/08/17): fix benthic submodel 
    6468      !!---------------------------------------------------------------------- 
    6569      !! Process benthic in/out fluxes 
    6670      !! These can be handled outside of the 3D calculations since the 
    67       !! benthic pools (and fluxes) are 2D in nature; this code is 
    68       !! (shamelessly) borrowed from corresponding code in the LOBSTER 
    69       !! model 
     71      !! benthic pools (and fluxes) are 2D in nature; this code was 
     72      !! developed with help from George Nurser (NOC); it cannot be run 
     73      !! in a configuration with variable time-stepping with depth 
    7074      !!---------------------------------------------------------------------- 
    7175      !! 
    72       !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt 
     76      !! time-step calculation 
    7377      if (jorgben.eq.1) then 
    74          za_sed_n(:,:)  = zn_sed_n(:,:)  +                                 &  
    75                           ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -          & 
    76                             f_benout_n(:,:)  ) * (rdt / 86400.) 
     78         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * & 
     79                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  )) 
     80         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * & 
     81                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) )) 
     82         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * & 
     83                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  )) 
     84      endif 
     85      if (jinorgben.eq.1) then 
     86         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * & 
     87                          ( f_fbenin_si(:,:) - f_benout_si(:,:) )) 
     88         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * & 
     89                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) )) 
     90      endif 
     91      !! 
     92      !! time-level calculation 
     93      if (jorgben.eq.1) then 
     94         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * & 
     95                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  )) 
    7796         zn_sed_n(:,:)  = za_sed_n(:,:) 
    78          !! 
    79          za_sed_fe(:,:) = zn_sed_fe(:,:) +                                 & 
    80                           ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -          & 
    81                             f_benout_fe(:,:) ) * (rdt / 86400.) 
     97         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * & 
     98                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) )) 
    8299         zn_sed_fe(:,:) = za_sed_fe(:,:) 
    83          !! 
    84          za_sed_c(:,:)  = zn_sed_c(:,:)  +                                 & 
    85                           ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -          & 
    86                             f_benout_c(:,:)  ) * (rdt / 86400.) 
     100         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * & 
     101                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  )) 
    87102         zn_sed_c(:,:)  = za_sed_c(:,:) 
    88103      endif 
    89104      if (jinorgben.eq.1) then 
    90          za_sed_si(:,:) = zn_sed_si(:,:) +                                 &  
    91                           ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        & 
    92                           (rdt / 86400.) 
     105         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * & 
     106                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) )) 
    93107         zn_sed_si(:,:) = za_sed_si(:,:) 
    94          !! 
    95          za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 & 
    96                           ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        & 
    97                           (rdt / 86400.) 
     108         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * & 
     109                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) )) 
    98110         zn_sed_ca(:,:) = za_sed_ca(:,:) 
    99111      endif 
    100       !! 
    101       if (ibenthic.eq.2) then 
    102          !! The code below (in this if ... then ... endif loop) is 
    103          !! effectively commented out because it does not work as  
    104          !! anticipated; it can be deleted at a later date 
    105          if (jorgben.eq.1) then 
    106             za_sed_n(:,:)  = ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  -       & 
    107                                f_benout_n(:,:)  ) * rdt 
    108             za_sed_fe(:,:) = ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) -       & 
    109                                f_benout_fe(:,:) ) * rdt 
    110             za_sed_c(:,:)  = ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  -       & 
    111                                f_benout_c(:,:)  ) * rdt 
    112          endif 
    113          if (jinorgben.eq.1) then 
    114             za_sed_si(:,:) = ( f_fbenin_si(:,:) - f_benout_si(:,:) ) * rdt 
    115             za_sed_ca(:,:) = ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) * rdt 
    116          endif 
    117          !! 
    118          !! Leap-frog scheme - only in explicit case, otherwise the  
    119          !! time stepping is already being done in trczdf 
    120          !! IF( l_trczdf_exp .AND. (ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN 
    121          !!    zfact = 2. * rdttra(jk) * FLOAT( ndttrc ) 
    122          !!    IF( neuler == 0 .AND. kt == nittrc000 )  zfact = rdttra(jk) *  
    123          !!                                             FLOAT(ndttrc) 
    124          !!    if (jorgben.eq.1) then 
    125          !!       za_sed_n(:,:)  = zb_sed_n(:,:)  + ( zfact * za_sed_n(:,:)  ) 
    126          !!      za_sed_fe(:,:) = zb_sed_fe(:,:) + ( zfact * za_sed_fe(:,:) ) 
    127          !!       za_sed_c(:,:)  = zb_sed_c(:,:)  + ( zfact * za_sed_c(:,:)  ) 
    128          !!    endif 
    129          !!    if (jinorgben.eq.1) then 
    130          !!       za_sed_si(:,:) = zb_sed_si(:,:) + ( zfact * za_sed_si(:,:) ) 
    131          !!       za_sed_ca(:,:) = zb_sed_ca(:,:) + ( zfact * za_sed_ca(:,:) ) 
    132          !!    endif 
    133          !! ENDIF 
    134          !!  
    135          !! Time filter and swap of arrays 
    136          IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd  ) THEN ! centred or tvd scheme 
    137             IF( neuler == 0 .AND. kt == nittrc000 ) THEN 
    138                if (jorgben.eq.1) then 
    139                   zb_sed_n(:,:)  = zn_sed_n(:,:) 
    140                   zn_sed_n(:,:)  = za_sed_n(:,:) 
    141                   za_sed_n(:,:)  = 0.0 
    142                   !! 
    143                   zb_sed_fe(:,:) = zn_sed_fe(:,:) 
    144                   zn_sed_fe(:,:) = za_sed_fe(:,:) 
    145                   za_sed_fe(:,:) = 0.0 
    146                   !! 
    147                   zb_sed_c(:,:)  = zn_sed_c(:,:) 
    148                   zn_sed_c(:,:)  = za_sed_c(:,:) 
    149                   za_sed_c(:,:)  = 0.0 
    150                endif 
    151                if (jinorgben.eq.1) then 
    152                   zb_sed_si(:,:) = zn_sed_si(:,:) 
    153                   zn_sed_si(:,:) = za_sed_si(:,:) 
    154                   za_sed_si(:,:) = 0.0 
    155                   !! 
    156                   zb_sed_ca(:,:) = zn_sed_ca(:,:) 
    157                   zn_sed_ca(:,:) = za_sed_ca(:,:) 
    158                   za_sed_ca(:,:) = 0.0 
    159                endif 
    160             ELSE 
    161                if (jorgben.eq.1) then 
    162                   zb_sed_n(:,:)  = (atfp  *                                 & 
    163                                     ( zb_sed_n(:,:)  + za_sed_n(:,:)  )) +  & 
    164                                       (atfp1 * zn_sed_n(:,:) ) 
    165                   zn_sed_n(:,:)  = za_sed_n(:,:) 
    166                   za_sed_n(:,:)  = 0.0 
    167                   !! 
    168                   zb_sed_fe(:,:) = (atfp  *                                 & 
    169                                     ( zb_sed_fe(:,:) + za_sed_fe(:,:) )) +  & 
    170                                       (atfp1 * zn_sed_fe(:,:)) 
    171                   zn_sed_fe(:,:) = za_sed_fe(:,:) 
    172                   za_sed_fe(:,:) = 0.0 
    173                   !! 
    174                   zb_sed_c(:,:)  = (atfp  *                                 & 
    175                                     ( zb_sed_c(:,:)  + za_sed_c(:,:)  )) +  & 
    176                                       (atfp1 * zn_sed_c(:,:) ) 
    177                   zn_sed_c(:,:)  = za_sed_c(:,:) 
    178                   za_sed_c(:,:)  = 0.0 
    179                endif 
    180                if (jinorgben.eq.1) then 
    181                   zb_sed_si(:,:) = (atfp  *                                 & 
    182                                     ( zb_sed_si(:,:) + za_sed_si(:,:) )) +  & 
    183                                       (atfp1 * zn_sed_si(:,:)) 
    184                   zn_sed_si(:,:) = za_sed_si(:,:) 
    185                   za_sed_si(:,:) = 0.0 
    186                   !! 
    187                   zb_sed_ca(:,:) = (atfp  *                                 & 
    188                                     ( zb_sed_ca(:,:) + za_sed_ca(:,:) )) +  & 
    189                                       (atfp1 * zn_sed_ca(:,:)) 
    190                   zn_sed_ca(:,:) = za_sed_ca(:,:) 
    191                   za_sed_ca(:,:) = 0.0 
    192                endif 
    193             ENDIF 
    194          ELSE                   !  case of smolar scheme or muscl 
    195             if (jorgben.eq.1) then 
    196                zb_sed_n(:,:)  = za_sed_n(:,:) 
    197                zn_sed_n(:,:)  = za_sed_n(:,:) 
    198                za_sed_n(:,:)  = 0.0 
    199                !! 
    200                zb_sed_fe(:,:) = za_sed_fe(:,:) 
    201                zn_sed_fe(:,:) = za_sed_fe(:,:) 
    202                za_sed_fe(:,:) = 0.0 
    203                !! 
    204                zb_sed_c(:,:)  = za_sed_c(:,:) 
    205                zn_sed_c(:,:)  = za_sed_c(:,:) 
    206                za_sed_c(:,:)  = 0.0 
    207             endif 
    208             if (jinorgben.eq.1) then 
    209                zb_sed_si(:,:) = za_sed_si(:,:) 
    210                zn_sed_si(:,:) = za_sed_si(:,:) 
    211                za_sed_si(:,:) = 0.0 
    212                !! 
    213                zb_sed_ca(:,:) = za_sed_ca(:,:) 
    214                zn_sed_ca(:,:) = za_sed_ca(:,:) 
    215                za_sed_ca(:,:) = 0.0 
    216             endif 
    217          ENDIF 
    218       endif 
     112# endif       
    219113 
    220114#  if defined key_debug_medusa 
     
    253147                  fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) 
    254148                  fq2 = fq0 + fq1 
    255                   IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               & 
    256                      'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err)=',      & 
    257                      fq0,fq1,fq2 
     149                  fq3 = f_benout_n(ji,jj) 
     150                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     151                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      & 
     152                     fq0,fq1,fq2,fq3 
    258153               ENDIF 
    259154            ENDDO 
     
    266161                  fq1 = f_fbenin_si(ji,jj) 
    267162                  fq2 = fq0 + fq1 
    268                   IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               & 
    269                      'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err)=',     & 
    270                      fq0,fq1,fq2 
     163                  fq3 = f_benout_si(ji,jj) 
     164                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     165                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     166                     fq0,fq1,fq2,fq3 
    271167               ENDIF 
    272168            ENDDO 
     
    278174                  fq0 = fflx_c(ji,jj) 
    279175                  fq1 = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) + f_fbenin_ca(ji,jj) 
    280                   fq2 = f_co2flux(ji,jj) * e3t_n(ji,jj,1) 
     176                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1) 
    281177                  fq3 = fq0 + fq1 
    282                   IF (lwp) write (numout,'(a,2i3,a,4f15.10)')               & 
    283                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err)=',  & 
    284                     fq0,fq1,fq2,fq3 
    285                ENDIF 
    286             ENDDO 
    287          ENDDO    
    288          !! alkalinity 
    289          DO jj = 2,jpjm1 
    290             DO ji = 2,jpim1 
    291                if (tmask(ji,jj,1) == 1) then 
    292                   fq0 = fflx_a(ji,jj) 
    293                   fq1 = 2.0 * f_fbenin_ca(ji,jj) 
    294                   fq2 = fq0 + fq1 
    295                   IF (lwp) write (numout,'(a,2i3,a,3f15.10)')               & 
    296                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err)=',     & 
    297                      fq0,fq1,fq2 
     178                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj) 
     179                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   & 
     180                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', & 
     181                     fq0,fq1,fq2,fq3,fq4 
     182                ENDIF 
     183             ENDDO 
     184          ENDDO    
     185          !! alkalinity 
     186          DO jj = 2,jpjm1 
     187             DO ji = 2,jpim1 
     188                if (tmask(ji,jj,1) == 1) then 
     189                   fq0 = fflx_a(ji,jj) 
     190                   fq1 = 2.0 * f_fbenin_ca(ji,jj) 
     191                   fq2 = fq0 + fq1 
     192                   fq3 = 2.0 * f_benout_ca(ji,jj) 
     193                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     194                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     195                      fq0,fq1,fq2,fq3 
    298196               ENDIF 
    299197            ENDDO 
     
    360258            CALL iom_put( "OCAL_LVL"  , fccd ) 
    361259         ENDIF 
     260         IF ( med_diag%CHL_MLD%dgsave ) THEN 
     261            CALL iom_put( "CHL_MLD"  , fchl_ml ) 
     262         ENDIF 
    362263         IF ( med_diag%PN_JLIM%dgsave ) THEN 
    363264            CALL iom_put( "PN_JLIM"  , fjln2d ) 
Note: See TracChangeset for help on using the changeset viewer.