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 9157 for branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90 – NEMO

Ignore:
Timestamp:
2017-12-21T16:51:24+01:00 (6 years ago)
Author:
jpalmier
Message:

JPALM - merge GO6 branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_GO6_Carb_Fail_from_GO6_8356/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r9070 r9157  
    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 
     
    3839# endif 
    3940      USE lbclnk,            ONLY: lbc_lnk 
     41      USE oce,               ONLY: chloro_out_cpl  
    4042      USE par_medusa,        ONLY: jp_medusa_2d, jp_medusa_3d,          & 
    41                                    jp_medusa_trd 
     43                                   jp_medusa_trd, jpchd, jpchn 
    4244      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk 
    4345      USE phycst,            ONLY: rsmall 
     46      USE sbc_oce,           ONLY: lk_oasis 
    4447      USE sms_medusa,        ONLY: jinorgben, jorgben,                  & 
    4548                                   f3_co3, f3_h2co3, f3_hco3,           & 
     
    5053                                   zb_sed_n, zb_sed_si,                 & 
    5154                                   zn_sed_c, zn_sed_ca, zn_sed_fe,      & 
    52                                    zn_sed_n, zn_sed_si 
    53       USE trc,               ONLY: med_diag, nittrc000  
     55                                   zn_sed_n, zn_sed_si, zn_chl_srf,     & 
     56                                   scl_chl, chl_out 
     57      USE trc,               ONLY: med_diag, nittrc000, trn  
    5458      USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd 
    5559  
     
    6266      REAL(wp) :: fq0,fq1,fq2,fq3 
    6367 
     68# if defined key_roam                      
     69      !!---------------------------------------------------------------------- 
     70      !! AXY (09/08/17): fix benthic submodel 
    6471      !!---------------------------------------------------------------------- 
    6572      !! Process benthic in/out fluxes 
    6673      !! 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 
     74      !! benthic pools (and fluxes) are 2D in nature; this code was 
     75      !! developed with help from George Nurser (NOC); it cannot be run 
     76      !! in a configuration with variable time-stepping with depth 
    7077      !!---------------------------------------------------------------------- 
    7178      !! 
    72       !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt 
     79      !! time-step calculation 
    7380      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.) 
     81         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * & 
     82                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  )) 
     83         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * & 
     84                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) )) 
     85         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * & 
     86                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  )) 
     87      endif 
     88      if (jinorgben.eq.1) then 
     89         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * & 
     90                          ( f_fbenin_si(:,:) - f_benout_si(:,:) )) 
     91         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * & 
     92                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) )) 
     93      endif 
     94      !! 
     95      !! time-level calculation 
     96      if (jorgben.eq.1) then 
     97         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * & 
     98                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  )) 
    7799         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.) 
     100         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * & 
     101                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) )) 
    82102         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.) 
     103         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * & 
     104                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  )) 
    87105         zn_sed_c(:,:)  = za_sed_c(:,:) 
    88106      endif 
    89107      if (jinorgben.eq.1) then 
    90          za_sed_si(:,:) = zn_sed_si(:,:) +                                 &  
    91                           ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        & 
    92                           (rdt / 86400.) 
     108         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * & 
     109                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) )) 
    93110         zn_sed_si(:,:) = za_sed_si(:,:) 
    94          !! 
    95          za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 & 
    96                           ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        & 
    97                           (rdt / 86400.) 
     111         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * & 
     112                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) )) 
    98113         zn_sed_ca(:,:) = za_sed_ca(:,:) 
    99114      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 
     115# endif       
    219116 
    220117#  if defined key_debug_medusa 
     
    253150                  fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) 
    254151                  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 
     152                  fq3 = f_benout_n(ji,jj) 
     153                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     154                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      & 
     155                     fq0,fq1,fq2,fq3 
    258156               ENDIF 
    259157            ENDDO 
     
    266164                  fq1 = f_fbenin_si(ji,jj) 
    267165                  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 
     166                  fq3 = f_benout_si(ji,jj) 
     167                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     168                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     169                     fq0,fq1,fq2,fq3 
    271170               ENDIF 
    272171            ENDDO 
     
    278177                  fq0 = fflx_c(ji,jj) 
    279178                  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) 
     179                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1) 
    281180                  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 
     181                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj) 
     182                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   & 
     183                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', & 
     184                     fq0,fq1,fq2,fq3,fq4 
     185                ENDIF 
     186             ENDDO 
     187          ENDDO    
     188          !! alkalinity 
     189          DO jj = 2,jpjm1 
     190             DO ji = 2,jpim1 
     191                if (tmask(ji,jj,1) == 1) then 
     192                   fq0 = fflx_a(ji,jj) 
     193                   fq1 = 2.0 * f_fbenin_ca(ji,jj) 
     194                   fq2 = fq0 + fq1 
     195                   fq3 = 2.0 * f_benout_ca(ji,jj) 
     196                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     197                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     198                      fq0,fq1,fq2,fq3 
    298199               ENDIF 
    299200            ENDDO 
     
    333234            ENDDO 
    334235         ENDDO 
     236 
     237         !!!--------------------------------------------------------------- 
     238         !! Calculates Chl diag for UM coupling  
     239         !!!--------------------------------------------------------------- 
     240         !! JPALM -- 02-06-2017 -- 
     241         !! add Chl surf coupling 
     242         !! no need to output, just pass to cpl var 
     243         IF (lk_oasis) THEN 
     244            IF (chl_out.eq.1) THEN 
     245               !! export and scale surface chl 
     246               zn_chl_srf(:,:) = MAX( 0.0, (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6 ) 
     247                                 !! surf Chl in Kg-chl/m3 as needed for cpl 
     248            ELSEIF (chl_out.eq.2) THEN 
     249               !! export and scale mld chl 
     250               zn_chl_srf(:,:) = MAX( 0.0, fchl_ml(:,:) * 1.0E-6 ) 
     251                                 !! mld Chl in Kg-chl/m3 as needed for cpl 
     252            ENDIF 
     253            chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling Chl 
     254         END IF 
     255 
    335256         !!---------------------------------------------------------------- 
    336257         !! Add in XML diagnostics stuff 
     
    360281            CALL iom_put( "OCAL_LVL"  , fccd ) 
    361282         ENDIF 
     283         IF ( med_diag%CHL_MLD%dgsave ) THEN 
     284            CALL iom_put( "CHL_MLD"  , fchl_ml ) 
     285         ENDIF 
     286         IF (lk_oasis) THEN 
     287            IF ( med_diag%CHL_CPL%dgsave ) THEN 
     288               CALL iom_put( "CHL_CPL"  , chloro_out_cpl ) 
     289            ENDIF 
     290         ENDIF 
    362291         IF ( med_diag%PN_JLIM%dgsave ) THEN 
    363292            CALL iom_put( "PN_JLIM"  , fjln2d ) 
Note: See TracChangeset for help on using the changeset viewer.