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 10232 for branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

Ignore:
Timestamp:
2018-10-26T15:25:58+02:00 (6 years ago)
Author:
dford
Message:

Merge in revisions 8447:10159 of dev_r5518_GO6_package.

Location:
branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90

    r8442 r10232  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Add air-sea flux kill switch 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_medusa 
     
    6162# endif 
    6263                                   zchd, zchn, zdin, zsil 
    63       USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, tmask 
    64 # if defined key_iomput 
     64      USE dom_oce,           ONLY: e3t_0, gphit, tmask, mig, mjg 
     65# if defined key_vvl 
     66      USE dom_oce,           ONLY: e3t_n 
     67# endif 
    6568      USE iom,               ONLY: lk_iomput 
    66 # endif 
    6769      USE in_out_manager,    ONLY: lwp, numout 
    68       USE oce,               ONLY: PCO2a_in_cpl 
    6970      USE par_kind,          ONLY: wp 
    7071      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1 
    71       USE sbc_oce,           ONLY: fr_i, lk_oasis, qsr, wndm 
     72      USE sbc_oce,           ONLY: fr_i, qsr, wndm 
    7273      USE sms_medusa,        ONLY: jdms, jdms_input, jdms_model,          & 
    7374                                   jriver_alk, jriver_c,                  & 
    7475                                   jriver_n, jriver_si,                   & 
     76                                   ln_foam_medusa,                        & 
    7577                                   riv_alk, riv_c, riv_n, riv_si,         & 
    7678                                   zn_dms_chd, zn_dms_chn, zn_dms_din,    & 
    7779                                   zn_dms_mld, zn_dms_qsr,                & 
     80                                   f2_pco2w, f2_fco2w,                    & 
    7881                                   xnln, xnld  
    7982      USE trc,               ONLY: med_diag 
     
    8689#  else 
    8790      USE trcco2_medusa,     ONLY: trc_co2_medusa 
     91      USE mocsy_mainmod,     ONLY: p2fCO2 
    8892#  endif 
    8993      USE trcdms_medusa,     ONLY: trc_dms_medusa 
    9094      USE trcoxy_medusa,     ONLY: trc_oxy_medusa 
    9195# endif 
     96      USE lib_mpp,           ONLY: ctl_stop 
     97      USE trcstat,           ONLY: trc_rst_dia_stat  
    9298 
    9399   !!* Substitution 
     
    121127 
    122128# if defined key_roam 
     129      !! init 
     130      f_fco2w(:,:)       = 0.0 
     131      f_fco2atm(:,:)     = 0.0 
     132      f_schmidtco2(:,:)  = 0.0 
     133      f_kwco2(:,:)       = 0.0 
     134      f_co2starair(:,:)  = 0.0 
     135      f_dpco2(:,:)       = 0.0 
     136      f_rhosw(:,:)       = 0.0 
     137      f_K0(:,:)          = 0.0 
     138      !! air pressure (atm); ultimately this will use air  
     139      !! pressure at the base of the UKESM1 atmosphere  
     140      !!                                      
     141      f_pp0(:,:)   = 1.0 
     142 
     143 
    123144      !!----------------------------------------------------------- 
    124145      !! Air-sea gas exchange 
     
    133154         DO ji = 2,jpim1 
    134155            !! OPEN wet point IF..THEN loop 
    135             if (tmask(ji,jj,1) == 1) then 
    136                IF (lk_oasis) THEN 
    137                   !! use 2D atm xCO2 from atm coupling 
    138                   f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) 
    139                ENDIF 
     156            IF (tmask(ji,jj,1) == 1) then 
    140157               !! 
    141158               !! AXY (23/06/15): as part of an effort to update the  
     
    161178               'air-sea: carb-chem kt = ', kt 
    162179               CALL flush(numout) 
     180               !! JPALM add carb print: 
     181               call trc_rst_dia_stat(f_xco2a(:,:), 'f_xco2a') 
     182               call trc_rst_dia_stat(wndm(:,:), 'wndm') 
     183               call trc_rst_dia_stat(f_kw660(:,:), 'f_kw660') 
     184               call trc_rst_dia_stat(ztmp(:,:), 'ztmp') 
     185               call trc_rst_dia_stat(zsal(:,:), 'zsal') 
     186               call trc_rst_dia_stat(zalk(:,:), 'zalk') 
     187               call trc_rst_dia_stat(zdic(:,:), 'zdic') 
     188               call trc_rst_dia_stat(zsil(:,:), 'zsil') 
     189               call trc_rst_dia_stat(zpho(:,:), 'zpho') 
    163190#   endif 
     191#  if defined key_axy_carbchem 
     192#   if defined key_mocsy 
    164193      DO jj = 2,jpjm1 
    165194         DO ji = 2,jpim1 
    166195            if (tmask(ji,jj,1) == 1) then 
    167                !! air pressure (atm); ultimately this will use air  
    168                !! pressure at the base of the UKESM1 atmosphere  
    169                !!                                      
    170                f_pp0(ji,jj)   = 1.0 
    171                !! 
    172                !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp(ji,jj) 
    173                !! IF(lwp) WRITE(numout,*) ' MEDUSA wndm    =', wndm(ji,jj) 
    174                !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj) 
    175                !! 
    176 #  if defined key_axy_carbchem 
    177 #   if defined key_mocsy 
     196               !! 
     197               !! Jpalm -- 12-09-2017 -- add extra check after reccurent 
     198               !!          carbonate failure in the coupled run. 
     199               !!          must be associated to air-sea flux or air xCO2... 
     200               !!          Check MOCSY inputs 
     201               IF ( (zsal(ji,jj) > 75.0 ).OR.(zsal(ji,jj) < 0.0 ) .OR.        & 
     202                    (ztmp(ji,jj) > 50.0 ).OR.(ztmp(ji,jj) < -20.0 ) .OR.      & 
     203                    (zalk(ji,jj) > 35.0E2 ).OR.(zalk(ji,jj) <= 0.0 ) .OR.     & 
     204                    (zdic(ji,jj) > 35.0E2 ).OR.(zdic(ji,jj) <= 0.0 ) .OR.     & 
     205                    (f_kw660(ji,jj) > 1.0E-2 ).OR.(f_kw660(ji,jj) < 0.0 ) ) THEN 
     206                  IF(lwp) THEN  
     207                      WRITE(numout,*) ' surface T = ',ztmp(ji,jj) 
     208                      WRITE(numout,*) ' surface S = ',zsal(ji,jj) 
     209                      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj) 
     210                      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj) 
     211                      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj) 
     212                      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)    
     213                      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj) 
     214                      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj) 
     215                      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj) 
     216                      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj) 
     217                      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj) 
     218                      WRITE(numout,*) ' MOCSY input: ji =', mig(ji),' jj = ', mjg(jj),  & 
     219                                       ' kt = ', kt  
     220                      WRITE(numout,*) 'MEDUSA - Air-Sea INPUT: unrealistic surface Carb. Chemistry' 
     221                  ENDIF      
     222                  CALL ctl_stop( 'MEDUSA - Air-Sea INPUT: ',             & 
     223                                 'unrealistic surface Carb. Chemistry -- INPUTS' ) 
     224               ENDIF      
    178225               !! 
    179226               !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
     
    200247               f_TALK(ji,jj) = (zalk(ji,jj) / f_rhosw(ji,jj)) * 1000. 
    201248               f_dcf(ji,jj)  = f_rhosw(ji,jj) 
     249               !! Jpalm -- 12-09-2017 -- add extra check after reccurent 
     250               !!          carbonate failure in the coupled run. 
     251               !!          must be associated to air-sea flux or air xCO2... 
     252               !!          Check MOCSY outputs 
     253               !!=================== 
     254               !! Jpalm -- 19-02-2018 -- remove the cap - only check MOCSY inputs 
     255               !!       because of specific area in arabic sea where strangely 
     256               !!       with core 2 forcing, ALK is lower than DIC and result in  
     257               !!       Enormous dpco2 - even if all carb chem caract are OK. 
     258               !!       and this check stops the model. 
     259               !!       --Input checks are already more than enough to stop the 
     260               !!       model if carb chem goes crazy.  
     261               !!       we remove the mocsy output checks 
     262               !!=================== 
     263               !!IF ( (f_pco2w(ji,jj) > 1.E4 ).OR.(f_pco2w(ji,jj) < 0.0 ) .OR.     & 
     264               !!    (f_fco2w(ji,jj) > 1.E4 ).OR.(f_fco2w(ji,jj) < 0.0 ) .OR.     &    
     265               !!    (f_fco2atm(ji,jj) > 1.E4 ).OR.(f_fco2atm(ji,jj) < 0.0 ) .OR.     & 
     266               !!    (f_co2flux(ji,jj) > 1.E-1 ).OR.(f_co2flux(ji,jj) < -1.E-1 ) .OR.     & 
     267               !!    (f_dpco2(ji,jj) > 1.E4 ).OR.(f_dpco2(ji,jj) < -1.E4 ) ) THEN 
     268               !!  IF(lwp) THEN  
     269               !!      WRITE(numout,*) ' surface T = ',ztmp(ji,jj) 
     270               !!      WRITE(numout,*) ' surface S = ',zsal(ji,jj) 
     271               !!      WRITE(numout,*) ' surface ALK = ',zalk(ji,jj) 
     272               !!      WRITE(numout,*) ' surface DIC = ',zdic(ji,jj) 
     273               !!      WRITE(numout,*) ' KW660 = ',f_kw660(ji,jj) 
     274               !!      WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj)    
     275               !!      WRITE(numout,*) ' surface pco2w  = ',f_pco2w(ji,jj) 
     276               !!      WRITE(numout,*) ' surface fco2w  = ',f_fco2w(ji,jj) 
     277               !!      WRITE(numout,*) ' surface fco2a  = ',f_fco2atm(ji,jj) 
     278               !!      WRITE(numout,*) ' surface co2flx = ',f_co2flux(ji,jj) 
     279               !!      WRITE(numout,*) ' surface dpco2  = ',f_dpco2(ji,jj) 
     280               !!      WRITE(numout,*) ' MOCSY output: ji =', mig(ji),' jj = ', mjg(jj),  & 
     281               !!                        ' kt = ', kt      
     282               !!      WRITE(numout,*) 'MEDUSA - Air-Sea OUTPUT: unrealistic surface Carb. Chemistry' 
     283               !!  ENDIF      
     284               !!  CALL ctl_stop( 'MEDUSA - Air-Sea OUTPUT: ',            & 
     285               !!                 'unrealistic surface Carb. Chemistry -- OUTPUTS' ) 
     286               !!ENDIF      
    202287            ENDIF 
    203288         ENDDO 
    204289      ENDDO 
    205290 
     291#   if defined key_debug_medusa 
     292               !! JPALM add carb print: 
     293               call trc_rst_dia_stat(f_pco2w(:,:), 'f_pco2w') 
     294               call trc_rst_dia_stat(f_fco2w(:,:), 'f_fco2w') 
     295               call trc_rst_dia_stat(f_fco2atm(:,:), 'f_fco2atm') 
     296               call trc_rst_dia_stat(f_schmidtco2(:,:), 'f_schmidtco2') 
     297               call trc_rst_dia_stat(f_kwco2(:,:), 'f_kwco2') 
     298               call trc_rst_dia_stat(f_co2starair(:,:), 'f_co2starair') 
     299               call trc_rst_dia_stat(f_co2flux(:,:), 'f_co2flux') 
     300               call trc_rst_dia_stat(f_dpco2(:,:), 'f_dpco2') 
     301#   endif 
    206302#   else    
    207303 
     
    234330                     iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt 
    235331               endif 
     332               IF ( ln_foam_medusa ) THEN 
     333                  !! DAF (Aug 2017): calculate fCO2 for observation operator 
     334                  CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w ) 
     335               ENDIF 
    236336            ENDIF 
    237337         ENDDO 
     
    277377         ENDDO 
    278378      ENDDO 
     379#  endif 
     380 
     381#  if defined key_axy_killco2flux 
     382      !! AXY (18/08/17): single kill switch on air-sea CO2 flux for budget checking 
     383      f_co2flux(:,:) = 0. 
    279384#  endif 
    280385 
     
    408513                              CO2flux_conv 
    409514               !! ENDIF 
     515               IF ( ln_foam_medusa ) THEN 
     516                  !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator 
     517                  f2_pco2w(ji,jj) = f_pco2w(ji,jj) 
     518                  f2_fco2w(ji,jj) = f_pco2w(ji,jj) 
     519               ENDIF 
    410520               IF ( lk_iomput ) THEN 
    411521                  IF( med_diag%ATM_PCO2%dgsave ) THEN 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_med_diag_iomput.F90

    r8442 r10232  
    3131      !!------------------------------------------------------------------- 
    3232      USE bio_medusa_mod 
    33       USE dom_oce,           ONLY: e3t_0, e3t_n, mbathy, tmask 
     33      USE dom_oce,           ONLY: e3t_0, mbathy, tmask 
     34# if defined key_vvl 
     35      USE dom_oce,           ONLY: e3t_n 
     36# endif 
    3437      USE in_out_manager,    ONLY: lwp, numout 
    3538      USE par_oce,           ONLY: jpim1, jpjm1 
     
    718721CONTAINS 
    719722   SUBROUTINE bio_med_diag_iomput( )                    ! Empty routine 
     723      IMPLICIT NONE 
    720724      WRITE(*,*) 'bio_med_diag_iomput: You should not have seen this print! error?' 
    721725   END SUBROUTINE bio_med_diag_iomput 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag.F90

    r8442 r10232  
    3232      USE bio_med_diag_iomput_mod,  ONLY: bio_med_diag_iomput 
    3333      USE bio_medusa_mod 
    34       USE dom_oce,                  ONLY: e3t_0, e3t_n,                  & 
    35                                           gdepw_0, gdepw_n, tmask 
     34      USE dom_oce,                  ONLY: e3t_0, gdepw_0, tmask 
     35# if defined key_vvl 
     36      USE dom_oce,                  ONLY: e3t_n, gdepw_n 
     37# endif 
    3638      USE in_out_manager,           ONLY: lwp, numout 
    37 # if defined key_iomput 
    3839      USE iom,                      ONLY: lk_iomput 
    39 # endif 
    4040      USE par_oce,                  ONLY: jpim1, jpjm1 
    4141      USE sms_medusa,               ONLY: xrfn, xthetapd, xthetapn,      & 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_diag_slice.F90

    r8442 r10232  
    3535      USE dom_oce,           ONLY: tmask 
    3636      USE in_out_manager,    ONLY: lwp, numout 
    37 # if defined key_iomput 
    3837      USE iom,               ONLY: iom_put 
    39 # endif 
    4038      USE lbclnk,            ONLY: lbc_lnk 
    41       USE trc,               ONLY: trn 
    42       USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl 
     39      USE oce,               ONLY: CO2Flux_out_cpl, DMS_out_cpl 
    4340      USE par_oce,           ONLY: jpi, jpj 
    4441      USE sbc_oce,           ONLY: lk_oasis, qsr, wndm 
     
    4845                                   jdms, ocal_ccd, xpar, xze,              & 
    4946                                   zb_co2_flx, zb_dms_srf,                 & 
    50                                    zn_co2_flx, zn_dms_srf, zn_chl_srf 
     47                                   zn_co2_flx, zn_dms_srf 
    5148      USE trc,               ONLY: med_diag 
    5249 
     
    6663      !! 
    6764      IF (jk.eq.1) THEN 
    68          !! JPALM -- 02-06-2017 -- 
    69          !! add Chl surf coupling 
    70          !! no need to output, just pass to cpl var 
    71          IF (lk_oasis) THEN 
    72             zn_chl_srf(:,:) = (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6  !! surf Chl in Kg-chl/m3 as needed for cpl 
    73             chloro_out_cpl(:,:) = zn_chl_srf(:,:)        !! Coupling Chl 
    74          END IF 
    7565         IF( med_diag%MED_QSR%dgsave ) THEN 
    7666            CALL iom_put( "MED_QSR"  , qsr ) ! 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r8442 r10232  
    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 
     
    3233      !!---------------------------------------------------------------------- 
    3334      USE bio_medusa_mod 
    34       USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt, e3t_n, tmask 
     35      USE dom_oce,           ONLY: atfp, atfp1, neuler, rdt, tmask 
    3536      USE in_out_manager,    ONLY: lwp, numout 
    36 # if defined key_iomput 
    37       USE iom,               ONLY: iom_put, lk_iomput 
    38 # endif 
     37      USE iom,               ONLY: iom_put 
    3938      USE lbclnk,            ONLY: lbc_lnk 
     39      USE oce,               ONLY: chloro_out_cpl  
    4040      USE par_medusa,        ONLY: jp_medusa_2d, jp_medusa_3d,          & 
    41                                    jp_medusa_trd 
     41                                   jp_medusa_trd, jpchd, jpchn 
    4242      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk 
    4343      USE phycst,            ONLY: rsmall 
     44      USE sbc_oce,           ONLY: lk_oasis 
    4445      USE sms_medusa,        ONLY: jinorgben, jorgben,                  & 
    4546                                   f3_co3, f3_h2co3, f3_hco3,           & 
    4647                                   f3_omarg, f3_omcal, f3_pH,           & 
     48                                   ln_foam_medusa, mld_max, pgrow_avg,  & 
     49                                   ploss_avg, phyt_avg,                 & 
    4750                                   za_sed_c, za_sed_ca, za_sed_fe,      & 
    4851                                   za_sed_n, za_sed_si,                 & 
     
    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 
     59      USE zdfmxl,            ONLY: hmld 
    5560  
    5661      !! time (integer timestep) 
     
    6267      REAL(wp) :: fq0,fq1,fq2,fq3 
    6368 
     69# if defined key_roam                      
     70      !!---------------------------------------------------------------------- 
     71      !! AXY (09/08/17): fix benthic submodel 
    6472      !!---------------------------------------------------------------------- 
    6573      !! Process benthic in/out fluxes 
    6674      !! 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 
     75      !! benthic pools (and fluxes) are 2D in nature; this code was 
     76      !! developed with help from George Nurser (NOC); it cannot be run 
     77      !! in a configuration with variable time-stepping with depth 
    7078      !!---------------------------------------------------------------------- 
    7179      !! 
    72       !! IF(lwp) WRITE(numout,*) 'AXY: rdt = ', rdt 
     80      !! time-step calculation 
    7381      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.) 
     82         za_sed_n(:,:)  = zb_sed_n(:,:)  + ((2. * (rdt / 86400.)) * & 
     83                          ( f_sbenin_n(:,:)  + f_fbenin_n(:,:)  - f_benout_n(:,:)  )) 
     84         za_sed_fe(:,:) = zb_sed_fe(:,:) + ((2. * (rdt / 86400.)) * & 
     85                          ( f_sbenin_fe(:,:) + f_fbenin_fe(:,:) - f_benout_fe(:,:) )) 
     86         za_sed_c(:,:)  = zb_sed_c(:,:)  + ((2. * (rdt / 86400.)) * & 
     87                          ( f_sbenin_c(:,:)  + f_fbenin_c(:,:)  - f_benout_c(:,:)  )) 
     88      endif 
     89      if (jinorgben.eq.1) then 
     90         za_sed_si(:,:) = zb_sed_si(:,:) + ((2. * (rdt / 86400.)) * & 
     91                          ( f_fbenin_si(:,:) - f_benout_si(:,:) )) 
     92         za_sed_ca(:,:) = zb_sed_ca(:,:) + ((2. * (rdt / 86400.)) * & 
     93                          ( f_fbenin_ca(:,:) - f_benout_ca(:,:) )) 
     94      endif 
     95      !! 
     96      !! time-level calculation 
     97      if (jorgben.eq.1) then 
     98         zb_sed_n(:,:)  = zn_sed_n(:,:)  + (atfp * & 
     99                          ( za_sed_n(:,:)  - (2. * zn_sed_n(:,:))  + zb_sed_n(:,:)  )) 
    77100         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.) 
     101         zb_sed_fe(:,:) = zn_sed_fe(:,:) + (atfp * & 
     102                          ( za_sed_fe(:,:) - (2. * zn_sed_fe(:,:)) + zb_sed_fe(:,:) )) 
    82103         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.) 
     104         zb_sed_c(:,:)  = zn_sed_c(:,:)  + (atfp * & 
     105                          ( za_sed_c(:,:)  - (2. * zn_sed_c(:,:))  + zb_sed_c(:,:)  )) 
    87106         zn_sed_c(:,:)  = za_sed_c(:,:) 
    88107      endif 
    89108      if (jinorgben.eq.1) then 
    90          za_sed_si(:,:) = zn_sed_si(:,:) +                                 &  
    91                           ( f_fbenin_si(:,:) - f_benout_si(:,:) ) *        & 
    92                           (rdt / 86400.) 
     109         zb_sed_si(:,:) = zn_sed_si(:,:) + (atfp * & 
     110                          ( za_sed_si(:,:) - (2. * zn_sed_si(:,:)) + zb_sed_si(:,:) )) 
    93111         zn_sed_si(:,:) = za_sed_si(:,:) 
    94          !! 
    95          za_sed_ca(:,:) = zn_sed_ca(:,:) +                                 & 
    96                           ( f_fbenin_ca(:,:) - f_benout_ca(:,:) ) *        & 
    97                           (rdt / 86400.) 
     112         zb_sed_ca(:,:) = zn_sed_ca(:,:) + (atfp * & 
     113                          ( za_sed_ca(:,:) - (2. * zn_sed_ca(:,:)) + zb_sed_ca(:,:) )) 
    98114         zn_sed_ca(:,:) = za_sed_ca(:,:) 
    99115      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 
     116# endif       
     117 
     118      IF ( ln_foam_medusa ) THEN 
     119         !!---------------------------------------------------------------------- 
     120         !! Diagnostics required for ocean colour assimilation: 
     121         !! Mixed layer average phytoplankton growth, loss and concentration 
     122         !! Maximum mixed layer depth 
     123         !!---------------------------------------------------------------------- 
    117124         !! 
    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 
     125         DO jj = 2,jpjm1 
     126            DO ji = 2,jpim1 
     127               IF ( hmld(ji,jj) .GT. 0.0 ) THEN 
     128                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) / hmld(ji,jj) 
     129                  ploss_avg(ji,jj) = ploss_avg(ji,jj) / hmld(ji,jj) 
     130                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  / hmld(ji,jj) 
     131                  IF ( hmld(ji,jj) .GT. mld_max(ji,jj) ) THEN 
     132                     mld_max(ji,jj) = hmld(ji,jj) 
     133                  ENDIF 
     134               ENDIF 
     135            END DO 
     136         END DO 
     137      ENDIF 
    219138 
    220139#  if defined key_debug_medusa 
     
    253172                  fq1 = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) 
    254173                  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 
     174                  fq3 = f_benout_n(ji,jj) 
     175                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     176                     'AXY N   cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',      & 
     177                     fq0,fq1,fq2,fq3 
    258178               ENDIF 
    259179            ENDDO 
     
    266186                  fq1 = f_fbenin_si(ji,jj) 
    267187                  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 
     188                  fq3 = f_benout_si(ji,jj) 
     189                  if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     190                     'AXY Si  cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     191                     fq0,fq1,fq2,fq3 
    271192               ENDIF 
    272193            ENDDO 
     
    278199                  fq0 = fflx_c(ji,jj) 
    279200                  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) 
     201                  fq2 = f_co2flux(ji,jj) * fse3t(ji,jj,1) 
    281202                  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 
     203                  fq4 = f_benout_c(ji,jj) + f_benout_ca(ji,jj) 
     204                  if (lwp) write (numout,'a,2i3,a,5f15,5)')                   & 
     205                     'AXY C   cons: (i,j)=',ji,jj,', (flx,ben,asf,err,out)=', & 
     206                     fq0,fq1,fq2,fq3,fq4 
     207                ENDIF 
     208             ENDDO 
     209          ENDDO    
     210          !! alkalinity 
     211          DO jj = 2,jpjm1 
     212             DO ji = 2,jpim1 
     213                if (tmask(ji,jj,1) == 1) then 
     214                   fq0 = fflx_a(ji,jj) 
     215                   fq1 = 2.0 * f_fbenin_ca(ji,jj) 
     216                   fq2 = fq0 + fq1 
     217                   fq3 = 2.0 * f_benout_ca(ji,jj) 
     218                   if (lwp) write (numout,'a,2i3,a,4f15,5)')                   & 
     219                      'AXY alk cons: (i,j)=',ji,jj,', (flx,ben,err,out)=',     & 
     220                      fq0,fq1,fq2,fq3 
    298221               ENDIF 
    299222            ENDDO 
     
    333256            ENDDO 
    334257         ENDDO 
     258 
     259         !!!--------------------------------------------------------------- 
     260         !! Calculates Chl diag for UM coupling  
     261         !!!--------------------------------------------------------------- 
     262         !! JPALM -- 02-06-2017 -- 
     263         !! add Chl surf coupling 
     264         !! no need to output, just pass to cpl var 
     265         IF (lk_oasis) THEN 
     266            IF (chl_out.eq.1) THEN 
     267               !! export and scale surface chl 
     268               zn_chl_srf(:,:) = MAX( 0.0, (trn(:,:,1,jpchd) + trn(:,:,1,jpchn)) * 1.0E-6 ) 
     269                                 !! surf Chl in Kg-chl/m3 as needed for cpl 
     270            ELSEIF (chl_out.eq.2) THEN 
     271               !! export and scale mld chl 
     272               zn_chl_srf(:,:) = MAX( 0.0, fchl_ml(:,:) * 1.0E-6 ) 
     273                                 !! mld Chl in Kg-chl/m3 as needed for cpl 
     274            ENDIF 
     275            chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl        !! Coupling Chl 
     276         END IF 
     277 
    335278         !!---------------------------------------------------------------- 
    336279         !! Add in XML diagnostics stuff 
     
    360303            CALL iom_put( "OCAL_LVL"  , fccd ) 
    361304         ENDIF 
     305         IF ( med_diag%CHL_MLD%dgsave ) THEN 
     306            CALL iom_put( "CHL_MLD"  , fchl_ml ) 
     307         ENDIF 
     308         IF (lk_oasis) THEN 
     309            IF ( med_diag%CHL_CPL%dgsave ) THEN 
     310               CALL iom_put( "CHL_CPL"  , chloro_out_cpl ) 
     311            ENDIF 
     312         ENDIF 
    362313         IF ( med_diag%PN_JLIM%dgsave ) THEN 
    363314            CALL iom_put( "PN_JLIM"  , fjln2d ) 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90

    r8442 r10232  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Add slow-sinking detrius variables 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_medusa 
     
    3435      USE bio_medusa_mod 
    3536      USE par_oce,           ONLY: jpi, jpj, jpk 
    36       USE sms_medusa,        ONLY: jdms 
     37      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max 
    3738      USE trc,               ONLY: ln_diatrc, med_diag, nittrc000  
    38       USE in_out_manager,    ONLY: lwp 
    39  
    40 # if defined key_iomput 
    41       USE iom,               ONLY: lk_iomput, numout 
     39      USE in_out_manager,    ONLY: lwp, numout 
     40 
     41      USE iom,               ONLY: lk_iomput 
    4242      USE trcnam_medusa,     ONLY: trc_nam_iom_medusa 
    43 # endif 
    4443 
    4544      !! time (integer timestep) 
     
    161160      fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics 
    162161      fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics 
     162      !! AXY (16/08/17) 
     163      fchl_ml(:,:) = 0.0   !! mixed layer chlorophyll diagnostics 
    163164      !!  
    164165      fslownflux(:,:) = 0.0  
    165166      fslowcflux(:,:) = 0.0  
     167      !! 
     168      !! JPALM -- 21-09-2017 -- needed to debug air-sea carb 
     169      f_xco2a(:,:)  = 0.0 
     170      f_pco2w(:,:)  = 0.0 
     171      f_ph(:,:)     = 0.0 
     172      f_kw660(:,:)  = 0.0 
     173      ztmp(:,:)  = 0.0 
     174      zsal(:,:)  = 0.0 
     175      zalk(:,:)  = 0.0 
     176      zdic(:,:)  = 0.0 
     177      zsil(:,:)  = 0.0 
     178# if defined key_mocsy 
     179      ! zpho is only defined if key_mocsy 
     180      ! is active, so we must protect this 
     181      ! initialisation accordingly.  
     182      zpho(:,:)  = 0.0 
     183# endif 
     184      f_co2flux(:,:)  = 0.0  
     185      f_pco2atm(:,:)  = 0.0 
     186      f_h2co3(:,:)    = 0.0 
     187      f_hco3(:,:)     = 0.0 
     188      f_co3(:,:)      = 0.0 
     189      f_omarg(:,:)    = 0.0 
     190      f_omcal(:,:)    = 0.0 
     191      !! 
     192      !! AXY (08/08/17): zero slow detritus fluxes 
     193      fslowsink(:,:)  = 0.0 
     194# if defined key_roam 
     195      fslowsinkc(:,:) = 0.0 
     196# endif       
     197      !! 
     198      pgrow_avg(:,:) = 0.0 
     199      ploss_avg(:,:) = 0.0 
     200      phyt_avg(:,:)  = 0.0 
     201      IF( kt == nittrc000 ) THEN 
     202         mld_max(:,:) = 0.0 
     203      ENDIF 
    166204      !! 
    167205      !! allocate and initiate 2D diag 
     
    836874CONTAINS 
    837875   SUBROUTINE bio_medusa_init( )                   ! Empty routine 
     876      IMPLICIT NONE 
    838877      WRITE(*,*) 'bio_medusa_init: You should not have seen this print! error?' 
    839878   END SUBROUTINE bio_medusa_init 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_mod.F90

    r8442 r10232  
    77   !! History : 
    88   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     9   !!   -   ! 2017-08 (A. Yool)            Slow detritus, ML-avg chl variables 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_medusa 
     
    5455   !! AXY (01/03/10): add in mixed layer PP diagnostics 
    5556   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fprn_ml,fprd_ml 
     57   !! AXY (16/08/17): add in mixed layer chlorophyll diagnostic 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fchl_ml 
    5659   !! 
    5760   !! nutrient limiting factors 
     
    9497   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fregenfastc 
    9598# endif 
    96  
     99   !! 
     100   !! AXY (08/08/17): sinking of detritus moved here 
     101   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsink, fslowgain, fslowloss 
     102# if defined key_roam 
     103   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::    fslowsinkc, fslowgainc, fslowlossc 
     104# endif 
     105   !! 
    97106   !! Particle flux 
    98107   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: fdep1 
     
    287296               fjlim_pn(jpi,jpj),fjlim_pd(jpi,jpj),                   & 
    288297               fun_T(jpi,jpj),fun_Q10(jpi,jpj),                       & 
    289                fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),                     & 
     298               fprn_ml(jpi,jpj),fprd_ml(jpi,jpj),fchl_ml(jpi,jpj),    & 
    290299               fnln(jpi,jpj),ffln2(jpi,jpj),                          & 
    291300               fnld(jpi,jpj),ffld(jpi,jpj),fsld(jpi,jpj),             & 
     
    316325               fregenfastc(jpi,jpj),                                  & 
    317326# endif 
     327          fslowsink(jpi,jpj),fslowgain(jpi,jpj),                 & 
     328               fslowloss(jpi,jpj),                                    & 
     329# if defined key_roam 
     330          fslowsinkc(jpi,jpj),fslowgainc(jpi,jpj),               & 
     331               fslowlossc(jpi,jpj),                                   & 
     332# endif 
    318333               fdep1(jpi,jpj),                                        & 
    319334               ftempn(jpi,jpj),ftempsi(jpi,jpj),ftempfe(jpi,jpj),     & 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_update.F90

    r8442 r10232  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Amend slow-detritus bug 
     9   !!   -   ! 2017-08 (A. Yool)            Reformatting for clarity 
    810   !!---------------------------------------------------------------------- 
    911#if defined key_medusa 
     
    6062                                   fsil_cons, fsil_prod, fsdiss,             & 
    6163                                   ftempca, fthetad, fthetan,                & 
     64                                   fslowsink, fslowgain, fslowloss,          & ! AXY (22/08/17) 
     65                                   f_sbenin_n, f_sbenin_c,                   & 
    6266# if defined key_roam 
     67                                   fslowsinkc, fslowgainc, fslowlossc,       & ! AXY (22/08/17) 
    6368                                   fcar_cons, fcar_prod, fcomm_resp,         & 
    6469                                   fddc, fflx_a, fflx_c, fflx_o2, zoxy,      & 
     
    6671# endif 
    6772                                   zpds, zphd, zphn 
    68       USE dom_oce,           ONLY: e3t_0, e3t_n, gphit, mbathy, tmask 
     73      USE dom_oce,           ONLY: e3t_0, gphit, mbathy, tmask 
     74# if defined key_vvl 
     75      USE dom_oce,           ONLY: e3t_n 
     76# endif 
    6977      USE in_out_manager,    ONLY: lwp, numout 
    7078      USE lib_mpp,           ONLY: ctl_stop 
    7179      USE par_kind,          ONLY: wp 
    72       USE par_medusa,        ONLY: jp_medusa,                                & 
     80      USE par_medusa,        ONLY: jp_medusa, jp_msa0, jp_msa1,              & 
    7381                                   jpalk, jpchd, jpchn, jpdet, jpdic,        & 
    7482                                   jpdin, jpdtc, jpfer, jpoxy, jppds,        & 
     
    7886                                   jpoxy_lc, jppds_lc, jpphd_lc, jpphn_lc,   & 
    7987                                   jpsil_lc, jpzme_lc, jpzmi_lc 
    80       USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1 
     88      USE par_oce,           ONLY: jpi, jpim1, jpj, jpjm1, jpk 
    8189      USE par_trc,           ONLY: jptra 
    8290      USE sms_medusa,        ONLY: friver_dep,                               & 
     
    181189      ENDDO 
    182190 
    183       DO jj = 2,jpjm1 
    184          DO ji = 2,jpim1 
    185             if (tmask(ji,jj,jk) == 1) then 
    186                !! 
    187                !!---------------------------------------------------------- 
    188                !! detritus 
    189                btra(ji,jj,jpdet_lc) = b0 *                                      & 
    190                                           ! mort. losses  
    191                                    (fdpn(ji,jj) + ((1.0 - xfdfrac1) *        & 
    192                                                    fdpd(ji,jj)) +            & 
    193                                     fdzmi(ji,jj) +                           & 
    194                                     ((1.0 - xfdfrac2) * fdzme(ji,jj)) +      & 
    195                                           ! assim. inefficiency 
    196                                     ((1.0 - xbetan) * (finmi(ji,jj) +        & 
    197                                                        finme(ji,jj))) -      & 
    198                                           ! grazing and remin. 
    199                                     fgmid(ji,jj) - fgmed(ji,jj) -            & 
    200                                     fdd(ji,jj) +                             & 
    201                                           ! seafloor fast->slow 
    202                                     ffast2slown(ji,jj)) 
    203                !! 
     191      !!---------------------------------------------------------- 
     192      !! detritus 
     193      DO jj = 2,jpjm1 
     194         DO ji = 2,jpim1 
     195            if (tmask(ji,jj,jk) == 1) then 
     196               !! 
     197               btra(ji,jj,jpdet_lc) = b0 * (                           & 
     198                   fdpn(ji,jj)                                         & ! mort. losses  
     199                 + ((1.0 - xfdfrac1) * fdpd(ji,jj))                    & ! mort. losses  
     200                 + fdzmi(ji,jj)                                        & ! mort. losses 
     201                 + ((1.0 - xfdfrac2) * fdzme(ji,jj))                   & ! mort. losses 
     202                 + ((1.0 - xbetan) * (finmi(ji,jj) + finme(ji,jj)))    & ! assim. inefficiency 
     203                 - fgmid(ji,jj) - fgmed(ji,jj)                         & ! grazing 
     204                 - fdd(ji,jj)                                          & ! remin. 
     205                 + fslowgain(ji,jj) - fslowloss(ji,jj)                 & ! slow-sinking 
     206                 - (f_sbenin_n(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor 
     207                 + ffast2slown(ji,jj) )                                  ! seafloor fast->slow  
    204208            ENDIF 
    205209         ENDDO 
     
    266270                             ! mort. loss 
    267271                         ((1.0 - xfdfrac1) * fdpds(ji,jj)) +                 & 
    268                              ! egestion of grazed Si 
     272                             ! egestion of grazed Si 
    269273                         ((1.0 - xfdfrac3) * fgmepds(ji,jj)) +               & 
    270274                             ! fast diss. and metab. losses 
     
    305309                                          ffetop(ji,jj) + ffebot(ji,jj) -    & 
    306310                                          ffescav(ji,jj) ) 
     311            ENDIF 
     312         ENDDO 
     313      ENDDO 
     314 
    307315# if defined key_roam 
    308                !! 
    309                !!---------------------------------------------------------- 
    310                !! AXY (26/11/08): implicit detrital carbon change 
    311                btra(ji,jj,jpdtc_lc) = b0 * (                                    & 
    312                                             ! mort. losses 
    313                                          (xthetapn * fdpn(ji,jj)) +          & 
    314                                          ((1.0 - xfdfrac1) *                 & 
    315                                           (xthetapd * fdpd(ji,jj))) +        & 
    316                                          (xthetazmi * fdzmi(ji,jj)) +        & 
    317                                          ((1.0 - xfdfrac2) *                 & 
    318                                           (xthetazme * fdzme(ji,jj))) +      & 
    319                                              ! assim. inefficiency 
    320                                          ((1.0 - xbetac) *                   & 
    321                                           (ficmi(ji,jj) + ficme(ji,jj))) -   & 
    322                                              ! grazing and remin. 
    323                                          fgmidc(ji,jj) - fgmedc(ji,jj) -     & 
    324                                          fddc(ji,jj) +                       & 
    325                                              ! seafloor fast->slow 
    326                                          ffast2slowc(ji,jj) ) 
     316      !!---------------------------------------------------------- 
     317      !! AXY (26/11/08): implicit detrital carbon change 
     318      DO jj = 2,jpjm1 
     319         DO ji = 2,jpim1 
     320            if (tmask(ji,jj,jk) == 1) then  
     321               !! 
     322               btra(ji,jj,jpdtc_lc) = b0 * (                           & 
     323                   (xthetapn * fdpn(ji,jj))                            & ! mort. losses  
     324                 + ((1.0 - xfdfrac1) * (xthetapd * fdpd(ji,jj)))       & ! mort. losses  
     325                 + (xthetazmi * fdzmi(ji,jj))                          & ! mort. losses  
     326                 + ((1.0 - xfdfrac2) * (xthetazme * fdzme(ji,jj)))     & ! mort. losses  
     327                 + ((1.0 - xbetac) * (ficmi(ji,jj) + ficme(ji,jj)))    & ! assim. inefficiency 
     328                 - fgmidc(ji,jj) - fgmedc(ji,jj)                       & ! grazing 
     329                 - fddc(ji,jj)                                         & ! remin. 
     330                 + fslowgainc(ji,jj) - fslowlossc(ji,jj)               & ! slow-sinking 
     331                 - (f_sbenin_c(ji,jj) / fse3t(ji,jj,jk))               & ! slow-sinking loss to seafloor 
     332                 + ffast2slowc(ji,jj) )                                  ! seafloor fast->slow 
    327333            ENDIF 
    328334         ENDDO 
     
    575581                                                           f_o2flux(ji,jj)) 
    576582               endif 
     583            ENDIF 
     584         ENDDO 
     585      ENDDO 
    577586# endif 
    578             ENDIF 
    579          ENDDO 
    580       ENDDO 
    581587 
    582588# if defined key_debug_medusa 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/carb_chem.F90

    r8441 r10232  
    4040# endif 
    4141                                   zalk, zdic, zsal, zsil, ztmp  
    42       USE dom_oce,           ONLY: gdept_0, gdept_n, gdepw_0, gdepw_n,    & 
    43                                    gphit, mbathy, tmask 
     42      USE dom_oce,           ONLY: gdept_0, gdepw_0, gphit, mbathy, tmask 
     43# if defined key_vvl 
     44      USE dom_oce,           ONLY: gdept_n, gdepw_n 
     45# endif 
    4446      USE in_out_manager,    ONLY: lwp, numout 
    4547      USE oce,               ONLY: PCO2a_in_cpl, tsb, tsn 
     
    103105               !! OPEN wet point IF..THEN loop 
    104106               IF (tmask(ji,jj,jk).eq.1) THEN 
    105                   IF (lk_oasis) THEN 
    106                      !! use 2D atm xCO2 from atm coupling 
    107                      f_xco2a(ji,jj) = PCO2a_in_cpl(ji,jj) 
    108                   ENDIF 
    109107                  !! Do carbonate chemistry 
    110108                  !! 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus.F90

    r8441 r10232  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Revise slow-sinking of detritus 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_medusa 
     
    3536                                        f_sbenin_n, fdd,                   & 
    3637                                        idf, idfval,                       &    
    37 # if defined key_roam 
     38                                        fslowsink,                         & 
     39                                        fslowgain, fslowloss,              & 
     40# if defined key_roam 
     41                                        fslowsinkc,                        & 
     42                                        fslowgainc, fslowlossc,            & 
    3843                                        fddc,                              & 
    3944# endif 
    4045                                        fun_T, fun_Q10, zdet, zdtc 
    4146      USE detritus_fast_sink_mod, ONLY: detritus_fast_sink 
    42       USE dom_oce,                ONLY: mbathy, tmask 
     47      USE dom_oce,                ONLY: mbathy, e3t_0, gphit, tmask 
     48# if defined key_vvl 
     49      USE dom_oce,                ONLY: e3t_n 
     50# endif 
    4351      USE in_out_manager,         ONLY: lwp, numout 
    4452      USE par_oce,                ONLY: jpim1, jpjm1 
    4553      USE sms_medusa,             ONLY: jmd, jorgben, jsfd, vsed,          & 
    4654                                        xrfn, xmd, xmdc, xthetad 
     55 
     56   !!* Substitution 
     57#  include "domzgr_substitute.h90" 
    4758 
    4859      !! Level 
     
    123134         DO ji = 2,jpim1 
    124135            if (tmask(ji,jj,jk) == 1) then 
     136               !!---------------------------------------------------------------------- 
     137               !! Detritus sinking (AXY, 08/08/18) 
     138          !! Replaces slow-sinking done in trcsed_medusa.F90 
     139               !! 
     140               !! Uses the fslowsink variable to carry slow-sinking detritus from one 
     141               !! grid level to the next, variable fslowgain to "add" detritus sinking 
     142               !! from above and variable fslowloss to "subtract" detritus sinking out 
     143               !! to below; these variables appear in the differential equations of 
     144               !! detrital nitrogen and carbon below 
     145               !!---------------------------------------------------------------------- 
     146               !! 
     147               fslowgain(ji,jj)  = fslowsink(ji,jj) / fse3t(ji,jj,jk)                  ! = mmol N / m3 / d 
     148               if (jk.lt.mbathy(ji,jj)) then 
     149                  fslowloss(ji,jj)  = (zdet(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol N / m3 / d 
     150               else 
     151                  fslowloss(ji,jj)  = 0.                                               ! = mmol N / m3 / d 
     152               endif 
     153               fslowsink(ji,jj) = fslowloss(ji,jj) * fse3t(ji,jj,jk)                   ! = mmol N / m2 / d 
     154               !! 
     155#  if defined key_roam 
     156               fslowgainc(ji,jj) = fslowsinkc(ji,jj) / fse3t(ji,jj,jk)                 ! = mmol C / m3 / d 
     157               if (jk.lt.mbathy(ji,jj)) then 
     158                  fslowlossc(ji,jj) = (zdtc(ji,jj) * vsed * 86400.) / fse3t(ji,jj,jk)  ! = mmol C / m3 / d 
     159               else 
     160                  fslowlossc(ji,jj) = 0.                                               ! = mmol C / m3 / d 
     161               endif 
     162               fslowsinkc(ji,jj) = fslowlossc(ji,jj) * fse3t(ji,jj,jk)                 ! = mmol C / m2 / d 
     163#  endif 
     164            ENDIF 
     165         ENDDO 
     166      ENDDO 
     167 
     168      DO jj = 2,jpjm1 
     169         DO ji = 2,jpim1 
     170            if (tmask(ji,jj,jk) == 1) then 
    125171               !!--------------------------------------------------------- 
    126172               !! Detritus addition to benthos 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/detritus_fast_sink.F90

    r8441 r10232  
    6666                                   idf, idfval,                            & 
    6767                                   zdet, zdtc 
    68       USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,         & 
    69                                    gphit, mbathy, tmask 
     68      USE dom_oce,           ONLY: e3t_0, gdepw_0, gphit, mbathy, tmask 
     69# if defined key_vvl 
     70      USE dom_oce,           ONLY: e3t_n, gdepw_n 
     71# endif 
    7072      USE in_out_manager,    ONLY: lwp, numout 
    7173      USE oce,               ONLY: tsn 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/iron_chem_scav.F90

    r8441 r10232  
    3535                                   zdet, zfer, zphd, zphn, zzme, zzmi,    & 
    3636                                   idf, idfval                           
    37       USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n,        & 
    38                                    mbathy, tmask 
     37      USE dom_oce,           ONLY: e3t_0, gdepw_0, mbathy, tmask 
     38#if defined key_vvl 
     39      USE dom_oce,           ONLY: e3t_n, gdepw_n 
     40#endif 
    3941      USE par_kind,          ONLY: wp 
    4042      USE in_out_manager,    ONLY: lwp, numout 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90

    r8441 r10232  
    66   !! History : 
    77   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90 
     8   !!   -   ! 2017-08 (A. Yool)            Mean mixed layer chlorophyll 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_medusa 
     
    4243                                   zchd, zchn, zdet, zdin, zdtc,         & 
    4344                                   zfer, zpds, zphd, zphn, zsil,         & 
    44                                    zzme, zzmi 
    45       USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask 
     45                                   zzme, zzmi, fchl_ml 
     46      USE dom_oce,           ONLY: e3t_0, gdepw_0, tmask 
     47#if defined key_vvl 
     48      USE dom_oce,           ONLY: e3t_n, gdepw_n 
     49#endif 
    4650      USE in_out_manager,    ONLY: lwp, numout 
    4751      USE oce,               ONLY: tsn 
     
    5559                                   xvpd, xvpn, xxi 
    5660      USE zdfmxl,            ONLY: hmld 
     61      USE lbclnk,            ONLY: lbc_lnk 
    5762 
    5863   !!* Substitution 
     
    373378               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * & 
    374379                                                  fse3t(ji,jj,jk) * fq0) 
    375             ENDIF 
    376          ENDDO 
    377       ENDDO 
     380          !! AXY (16/08/17) 
     381          fchl_ml(ji,jj) = fchl_ml(ji,jj) + ((zchn(ji,jj) + zchd(ji,jj)) * & 
     382                                             (fse3t(ji,jj,jk) * fq0) / hmld(ji,jj)) 
     383            ENDIF 
     384         ENDDO 
     385      ENDDO 
     386      CALL lbc_lnk(fchl_ml(:,:),'T',1. ) 
    378387 
    379388      DO jj = 2,jpjm1 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90

    r8441 r10232  
    3636                                   fdpn, fdpn2, fdzme, fdzme2,             & 
    3737                                   fdzmi, fdzmi2, fsdiss, fsin,            & 
     38                                   fdep1, fprn, fprd,                      & 
     39                                   fgmepd, fgmepn, fgmipn,                 & 
    3840                                   zphd, zphn, zpds, zzme, zzmi 
    39       USE dom_oce,           ONLY: tmask 
     41      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask 
     42      USE par_kind,          ONLY: wp 
    4043      USE par_oce,           ONLY: jpim1, jpjm1 
    4144      USE phytoplankton_mod, ONLY: phytoplankton 
    4245      USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               & 
     46                                   ln_foam_medusa,                         & 
     47                                   pgrow_avg, ploss_avg, phyt_avg,         & 
    4348                                   xkphd, xkphn, xkzme, xkzmi,             & 
    4449                                   xmetapd, xmetapn, xmetazme, xmetazmi,   & 
    4550                                   xmpd, xmpn, xmzme, xmzmi, xsdiss 
     51      USE zdfmxl,            ONLY: hmld 
    4652      USE zooplankton_mod,   ONLY: zooplankton 
     53 
     54   !!* Substitution 
     55#  include "domzgr_substitute.h90" 
    4756 
    4857      !! Level 
     
    5059 
    5160      INTEGER :: ji, jj 
     61 
     62      REAL(wp) :: fq0 
    5263 
    5364      !!------------------------------------------------------------------- 
     
    188199      ENDDO 
    189200 
     201      IF ( ln_foam_medusa ) THEN 
     202         !! Mixed layer averages for ocean colour assimilation 
     203         !! 
     204         DO jj = 2,jpjm1 
     205            DO ji = 2,jpim1 
     206               IF (tmask(ji,jj,jk) == 1) THEN 
     207                  if (fdep1(ji,jj).le.hmld(ji,jj)) then 
     208                     !! this level is entirely in the mixed layer 
     209                     fq0 = 1.0 
     210                  elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then 
     211                     !! this level is entirely below the mixed layer 
     212                     fq0 = 0.0 
     213                  else 
     214                     !! this level straddles the mixed layer 
     215                     fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk) 
     216                  endif 
     217                  !! 
     218                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) +                   & 
     219                                     (((fprn(ji,jj) * zphn(ji,jj)) +      & 
     220                                       (fprd(ji,jj) * zphd(ji,jj))  ) *   & 
     221                                      fse3t(ji,jj,jk) * fq0) 
     222                  ploss_avg(ji,jj) = ploss_avg(ji,jj) +                   & 
     223                                     ((fgmepd(ji,jj) + fdpd(ji,jj) +      & 
     224                                       fdpd2(ji,jj)                +      & 
     225                                       fgmepn(ji,jj) + fdpn(ji,jj) +      & 
     226                                       fdpn2(ji,jj)  + fgmipn(ji,jj) ) *  & 
     227                                      fse3t(ji,jj,jk) * fq0) 
     228                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  +                   & 
     229                                     ((zphn(ji,jj) + zphd(ji,jj)) *       & 
     230                                      fse3t(ji,jj,jk) * fq0) 
     231               ENDIF 
     232            ENDDO 
     233         ENDDO 
     234      ENDIF 
     235 
    190236   END SUBROUTINE plankton 
    191237 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r8132 r10232  
    177177   INTEGER  ::  jdms_model   !: choice of DMS model passed to atmosphere 
    178178!!                              1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL 
    179 !! 
     179!! JPALM --19-12-2017 -- UM people need to tune the Anderson DMS scheme     
     180   REAL(wp) ::  dmsmin       !: Anderson DMS scheme - DMS minimum value 
     181   REAL(wp) ::  dmscut       !: Anderson DMS scheme - DMS cutoff value  
     182   REAL(wp) ::  dmsslp       !: Anderson DMS scheme - DMS slope  
     183!! FOR UKESM    
     184   REAL(wp) ::  scl_chl      !: scaling factor for tuned Chl passed to the UM  
     185   INTEGER  ::  chl_out      !: select Chl field exported and scaled for UM: 
     186                             !: 1- Surface Chl ; 2- MLD Chl 
    180187!! 
    181188   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   remdmp   !: depth dependent damping coefficient of passive tracers  
     
    205212   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth 
    206213!! 
     214!! 2D fields of pCO2 and fCO2 for observation operator 
     215   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_pco2w    !: 2D pCO2 
     216   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_fco2w    !: 2D fCO2 
     217!! 
    207218!! 2D fields of organic and inorganic material sedimented on the seafloor 
    208219   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before) 
     
    276287#if defined key_roam 
    277288!!---------------------------------------------------------------------- 
    278 !! Atmospheric pCO2 data (1859 to 2100 inclusive) 
    279 !!---------------------------------------------------------------------- 
    280 !! 
    281    REAL(wp), DIMENSION(242)         ::   hist_pco2 !: pCO2 
    282  
    283 # if defined key_rcp26 
    284       !! UKMO, run AJKKH + KAAEC, RCP 2.6, pCO2 time evolution 
    285       DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, & 
    286       & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, & 
    287       & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, & 
    288       & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, & 
    289       & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, & 
    290       & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, & 
    291       & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, & 
    292       & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, & 
    293       & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, & 
    294       & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, & 
    295       & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, & 
    296       & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, & 
    297       & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, & 
    298       & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, & 
    299       & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, & 
    300       & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, & 
    301       & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, & 
    302       & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, & 
    303       & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, & 
    304       & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, & 
    305       & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, & 
    306       & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, & 
    307       & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, & 
    308       & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, & 
    309       & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, & 
    310       & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, & 
    311       & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, & 
    312       & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, & 
    313       & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, & 
    314       & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, & 
    315       & 386.9310, 389.2150, 391.4910, 393.7710, 396.0460, & 
    316       & 398.3240, 400.6080, 402.8950, 405.1780, 407.4550, & 
    317       & 409.7260, 411.9930, 414.2500, 416.4410, 418.5280, & 
    318       & 420.5250, 422.4390, 424.2720, 426.0200, 427.6750, & 
    319       & 429.2360, 430.7050, 432.0850, 433.3580, 434.5140, & 
    320       & 435.5740, 436.5490, 437.4420, 438.2550, 438.9810, & 
    321       & 439.6110, 440.1430, 440.5770, 440.9450, 441.2660, & 
    322       & 441.5410, 441.7840, 442.0050, 442.2040, 442.3780, & 
    323       & 442.5210, 442.6200, 442.6720, 442.6810, 442.6540, & 
    324       & 442.5830, 442.4670, 442.3270, 442.1680, 441.9960, & 
    325       & 441.8060, 441.5930, 441.3440, 441.0540, 440.7230, & 
    326       & 440.3510, 439.9300, 439.4650, 438.9730, 438.4630, & 
    327       & 437.9400, 437.4020, 436.8400, 436.2640, 435.6850, & 
    328       & 435.1030, 434.5160, 433.9170, 433.3060, 432.7010, & 
    329       & 432.1110, 431.5380, 430.9810, 430.4320, 429.8860, & 
    330       & 429.3370, 428.7810, 428.2220, 427.6490, 427.0660, & 
    331       & 426.4890, 425.9270, 425.3840, 424.8610, 424.3540, & 
    332       & 423.8540, 423.3540, 422.8530, 422.3510, 421.8410, & 
    333       & 421.3250, 420.8190 / 
    334 # else 
    335       !! UKMO, run AJKKH + KAAEF, RCP 8.5, pCO2 time evolution 
    336       DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, & 
    337       & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, & 
    338       & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, & 
    339       & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, & 
    340       & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, & 
    341       & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, & 
    342       & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, & 
    343       & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, & 
    344       & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, & 
    345       & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, & 
    346       & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, & 
    347       & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, & 
    348       & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, & 
    349       & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, & 
    350       & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, & 
    351       & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, & 
    352       & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, & 
    353       & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, & 
    354       & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, & 
    355       & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, & 
    356       & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, & 
    357       & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, & 
    358       & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, & 
    359       & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, & 
    360       & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, & 
    361       & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, & 
    362       & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, & 
    363       & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, & 
    364       & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, & 
    365       & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, & 
    366       & 386.9420, 389.2540, 391.5670, 393.9370, 396.3920, & 
    367       & 398.9320, 401.5550, 404.2550, 407.0220, 409.8530, & 
    368       & 412.7470, 415.7050, 418.7210, 421.7880, 424.9180, & 
    369       & 428.1200, 431.3970, 434.7470, 438.1650, 441.6410, & 
    370       & 445.1700, 448.7530, 452.3920, 456.0950, 459.8810, & 
    371       & 463.7680, 467.7660, 471.8750, 476.0960, 480.4210, & 
    372       & 484.8390, 489.3470, 493.9430, 498.6400, 503.4380, & 
    373       & 508.3410, 513.3630, 518.5160, 523.8050, 529.2290, & 
    374       & 534.7780, 540.4450, 546.2230, 552.1120, 558.1110, & 
    375       & 564.2110, 570.4130, 576.7390, 583.1990, 589.7980, & 
    376       & 596.5390, 603.4110, 610.4060, 617.4940, 624.6500, & 
    377       & 631.8800, 639.1750, 646.5360, 653.9800, 661.5230, & 
    378       & 669.1840, 676.9570, 684.8290, 692.7790, 700.7690, & 
    379       & 708.8050, 716.8870, 725.0020, 733.1770, 741.3900, & 
    380       & 749.6700, 758.0480, 766.5050, 775.0350, 783.6110, & 
    381       & 792.2200, 800.8740, 809.5680, 818.2760, 827.0090, & 
    382       & 835.8020, 844.6550, 853.5730, 862.5690, 871.6190, & 
    383       & 880.7020, 889.8240, 898.9590, 908.1270, 917.3080, & 
    384       & 926.4960, 935.7040 / 
    385 # endif 
     289!! JPALM -- change hist_pco2 init 
     290!!---------------------------------------------------------------------- 
     291!! 
     292   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   hist_pco2 !: pCO2 
     293   INTEGER  :: co2_rec 
     294   REAL(wp) :: co2_yinit, co2_yend    !: First and Last year read in the xCO2.atm file 
     295   REAL(wp) :: xobs_xco2a             !: Observed atmospheric xCO2, from namelist 
     296#endif 
     297 
     298!!---------------------------------------------------------------------- 
     299!! JPALM -- PI CO2 key 
     300!!---------------------------------------------------------------------- 
     301!! 
     302#if defined key_axy_pi_co2 
     303   LOGICAL , PUBLIC ::   lk_pi_co2 = .TRUE.  !: PI xCO2 used 
     304#else 
     305   LOGICAL , PUBLIC ::   lk_pi_co2 = .FALSE. !: PI xCO2 unused 
    386306#endif 
    387307 
     
    434354   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ??? 
    435355 
     356!!---------------------------------------------------------------------- 
     357!! Parameters required for ocean colour assimilation 
     358!!---------------------------------------------------------------------- 
     359!! 
     360   LOGICAL :: ln_foam_medusa                                 !: Flag to calculate and save diagnostics 
     361   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pgrow_avg  !: Mixed layer average phytoplankton growth 
     362   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss 
     363   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton 
     364   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth 
     365!! 
     366 
    436367   !!---------------------------------------------------------------------- 
    437368   !! NEMO/TOP 3.3 , NEMO Consortium (2010) 
     
    446377      !!---------------------------------------------------------------------- 
    447378      USE lib_mpp , ONLY: ctl_warn 
    448       INTEGER ::   ierr(8)        ! Local variables 
     379      INTEGER ::   ierr(9)        ! Local variables 
    449380      !!---------------------------------------------------------------------- 
    450381      ierr(:) = 0 
     
    456387      !* 2D and 3D fields of carbonate system parameters 
    457388      ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       & 
     389                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       & 
    458390         &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       & 
    459391         &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       & 
     
    504436         &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       & 
    505437         &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) )  
     438      !* Fields for ocean colour data assimilation 
     439      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       & 
     440         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) ) 
    506441#endif 
    507442      ! 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r8442 r10232  
    2323   !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6 
    2424   !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation 
     25   !!  -   !  2017-11  (J. Palm, A. Yool)     Diagnose tracer excursions 
    2526   !!---------------------------------------------------------------------- 
    2627   !! 
     
    7778                                            zfer, zpds, zphd, zphn, zsil,   & 
    7879                                            zzme, zzmi 
    79       USE dom_oce,                    ONLY: e3t_0, e3t_n,                   & 
    80                                             gdept_0, gdept_n,               & 
    81                                             gdepw_0, gdepw_n,               & 
    82                                             nday_year, nsec_day, nyear,     & 
    83                                             rdt, tmask 
     80      USE dom_oce,                    ONLY: e3t_0, gdept_0,  gdepw_0,       & 
     81                                            nday_year, nsec_day,            &  
     82                                            nyear, nyear_len, ndastp,       & 
     83                                            nsec_month,                     & 
     84                                            rdt, tmask, mig, mjg, nimpp,    & 
     85                                            njmpp  
     86#if defined key_vvl 
     87      USE dom_oce,                    ONLY: e3t_n, gdept_n,  gdepw_n 
     88#endif 
     89 
    8490      USE in_out_manager,             ONLY: lwp, numout, nn_date0 
    85 # if defined key_iomput 
    8691      USE iom,                        ONLY: lk_iomput 
    87 # endif 
    8892      USE lbclnk,                     ONLY: lbc_lnk 
    89       USE lib_mpp,                    ONLY: ctl_stop 
    90       USE oce,                        ONLY: tsb, tsn 
     93      USE lib_mpp,                    ONLY: mpp_max, mpp_maxloc,            &  
     94                                            mpp_min, mpp_minloc,            & 
     95                                            ctl_stop, ctl_warn, lk_mpp   
     96      USE oce,                        ONLY: tsb, tsn, PCO2a_in_cpl 
    9197      USE par_kind,                   ONLY: wp 
    9298      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     & 
     
    98104      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    99105      USE sbc_oce,                    ONLY: lk_oasis 
    100       USE sms_medusa,                 ONLY: hist_pco2 
     106      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend, & 
     107# if defined key_roam 
     108                                            xobs_xco2a,                     & 
     109# endif 
     110                                            pgrow_avg,                      & 
     111                                            ploss_avg, phyt_avg, mld_max,   & 
     112                                            lk_pi_co2, ln_foam_medusa 
    101113      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
    102114      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     
    110122      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice 
    111123      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin 
     124      USE trcstat,                    ONLY: trc_rst_dia_stat 
    112125 
    113126      IMPLICIT NONE 
     
    115128       
    116129      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90 
     130      PUBLIC   trc_bio_exceptionnal_fix ! here  
    117131 
    118132   !!* Substitution 
     
    176190      !!  
    177191      !! temporary variables 
    178       REAL(wp) ::    fq0,fq1,fq2,fq3,fq4 
     192      REAL(wp) ::    fq3,fq4 
     193      REAL(wp) ::    this_y, this_d, this_s, fyear 
     194      !! 
     195      !! T and S check temporary variable 
     196      REAL(wp) ::    sumtsn, tsnavg 
     197      INTEGER  ::    summask 
     198      CHARACTER(40) :: charout, charout2, charout3, charout4, charout5 
    179199      !! 
    180200      !!------------------------------------------------------------------ 
     
    283303      !!------------------------------------------------------------------ 
    284304      !! 
    285       !! what's atmospheric pCO2 doing? (data start in 1859) 
    286       iyr1  = nyear - 1859 + 1 
    287       iyr2  = iyr1 + 1 
    288       if (iyr1 .le. 1) then 
    289          !! before 1860 
    290          f_xco2a(:,:) = hist_pco2(1) 
    291       elseif (iyr2 .ge. 242) then 
    292          !! after 2099 
    293          f_xco2a(:,:) = hist_pco2(242) 
    294       else 
    295          !! just right 
    296          fq0 = hist_pco2(iyr1) 
    297          fq1 = hist_pco2(iyr2) 
    298          fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) 
    299          !! AXY (14/06/12): tweaked to make more sense (and be correct) 
    300 #  if defined key_bs_axy_yrlen 
    301          !! bugfix: for 360d year with HadGEM2-ES forcing 
    302          fq3 = (real(nday_year) - 1.0 + fq2) / 360.0   
    303 #  else 
    304          !! original use of 365 days (not accounting for leap year or  
    305          !! 360d year) 
    306          fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 
    307 #  endif 
    308          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
    309          f_xco2a(:,:) = fq4 
    310       endif 
    311 #  if defined key_axy_pi_co2 
    312       !! OCMIP pre-industrial pCO2 
    313       !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2 
    314       f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2  
    315 #  endif 
    316       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    317       !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day) 
    318       !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year) 
    319       !! AXY (29/01/14): remove surplus diagnostics 
    320       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0 
    321       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1 
    322       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2 
    323       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3 
    324       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1) 
     305      IF (lk_oasis) THEN  
     306         !! xCO2 from coupled 
     307         IF ( ( kt == nittrc000 ) .AND. lwp )     & 
     308             WRITE(numout,*) '** MEDUSA Atm xCO2 given by the UM **' 
     309         f_xco2a(:,:) = PCO2a_in_cpl(:,:) 
     310         !! Check the xCO2 from the UM is OK 
     311         !! piece of code moved from air-sea.F90 
     312         !!--- 
     313         DO jj = 2,jpjm1 
     314            DO ji = 2,jpim1 
     315               !! OPEN wet point IF..THEN loop 
     316               IF (tmask(ji,jj,1) == 1) then 
     317                  !!! Jpalm test on atm xCO2 
     318                  IF ( (f_xco2a(ji,jj) .GT. 10000.0 ).OR.     & 
     319                       (f_xco2a(ji,jj) .LE. 0.0 ) ) THEN 
     320                     IF(lwp) THEN 
     321                        WRITE(numout,*) ' atm xCO2 = ',f_xco2a(ji,jj),       & 
     322                                        ' -- ji =', mig(ji),' jj = ', mjg(jj) 
     323                     ENDIF 
     324                     CALL ctl_stop( 'MEDUSA - trc_bio :',     & 
     325                                    'unrealistic coupled atm xCO2 ' ) 
     326                  ENDIF 
     327               ENDIF 
     328            ENDDO 
     329         ENDDO 
     330         !!---  
     331      ELSEIF (lk_pi_co2) THEN 
     332         !! OCMIP pre-industrial xCO2 
     333         IF ( ( kt == nittrc000 ) .AND. lwp )     & 
     334             WRITE(numout,*) '** MEDUSA Atm xCO2 fixed to pre-industrial value **' 
     335         !! f_xco2a(:,:) = 284.725       !! CMIP5 pre-industrial pCO2 
     336         f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2  
     337      ELSEIF ( xobs_xco2a > 0.0 ) THEN 
     338         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a 
     339         f_xco2a(:,:) = xobs_xco2a 
     340      ELSE       
     341         !! xCO2 from file 
     342         !! AXY - JPALM new interpolation scheme usinf nyear_len 
     343         this_y = real(nyear) 
     344         this_d = real(nday_year) 
     345         this_s = real(nsec_day) 
     346         !! 
     347         fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1)) 
     348         !! 
     349         IF ( ( kt == nittrc000 ) .AND. lwp ) THEN  
     350            WRITE(numout,*) '** MEDUSA Atm xCO2 from file **' 
     351            WRITE(numout,*) ' MEDUSA year      =', this_y 
     352            WRITE(numout,*) ' Year length      =', real(nyear_len(1)) 
     353            WRITE(numout,*) ' MEDUSA nday_year =', this_d 
     354            WRITE(numout,*) ' MEDUSA nsec_day  =', this_s 
     355         ENDIF 
     356         !!  
     357         !! different case test  
     358         IF (fyear .LE. co2_yinit) THEN 
     359            !! before first record -- pre-industrial value 
     360            f_xco2a(:,:) = hist_pco2(1)   
     361         ELSEIF (fyear .GE. co2_yend) THEN 
     362            !! after last record - continue to use the last value 
     363            f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) ) 
     364         ELSE 
     365            !! just right 
     366            iyr1 = int(fyear - co2_yinit) + 1 
     367            iyr2 = iyr1 + 1  
     368            fq3 = fyear - real(iyr1) - co2_yinit + 1. 
     369            fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2)) 
     370            f_xco2a(:,:) = fq4 
     371            !! 
     372            IF ( ( kt == nittrc000 ) .AND. lwp ) THEN  
     373               WRITE(numout,*) ' MEDUSA year1     =', iyr1 
     374               WRITE(numout,*) ' MEDUSA year2     =', iyr2 
     375               WRITE(numout,*) ' xCO2 year1       =', hist_pco2(iyr1) 
     376               WRITE(numout,*) ' xCO2 year2       =', hist_pco2(iyr2) 
     377               WRITE(numout,*) ' Year2 weight     =', fq3 
     378            ENDIF 
     379         ENDIF 
     380      ENDIF 
     381  
     382      !! Writing xCO2 in output on start and on the 1st tsp of each  month 
     383      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     384           ( nsec_month .LE. INT(rdt) ) )  THEN 
     385           IF ( lwp )  WRITE(numout,*) ' *** Atm xCO2 *** -- kt:', kt,       & 
     386                                       ';  current date:', ndastp 
     387          call trc_rst_dia_stat(f_xco2a(:,:), 'atm xCO2') 
     388      ENDIF 
    325389# endif 
    326390 
     
    348412      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    349413      !!          ++ need to pass carb-chem output var through restarts 
    350       If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
    351            ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
     414      !!If ( (kt == nitt8rc000 .AND. .NOT.ln_rsttr) .OR.          & 
     415      !!     ( (mod(kt*rdt,2592000.)) == rdt) THEN 
     416      !!============================= 
     417      !! (Jpalm -- updated for restartability issues) 
     418      !! We want this to be start of month or if starting afresh from   
     419      !! climatology - marc 20/6/17  
     420      !!If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                         &  
     421      !!     ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN  
     422      !!============================= 
     423      !! Jpalm -- 15-02-2018 -- need to change 3D carb-chem call freq again. 
     424      !! previous call did not work, probably the (86400*mod(nn_date0,100) part 
     425      !! should not be in... 
     426      !! now use the NEMO calendar tool : nsec_month to be sure to call  
     427      !! at the beginning of a new month . 
     428      !! DAF: For FOAM we want to run daily 
     429      IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
     430           ( nsec_month .LE. INT(rdt) )          .OR.                        & 
     431           ( nsec_day   .LE. INT(rdt) .AND. ln_foam_medusa ) )  THEN 
     432           IF ( lwp )  WRITE(numout,*)                                       &          
     433                              ' *** 3D carb chem call *** -- kt:', kt,       & 
     434                              ';  current date:', ndastp  
    352435         !!--------------------------------------------------------------- 
    353436         !! Calculate the carbonate chemistry for the whole ocean on the first 
     
    450533 
    451534# if defined key_roam 
     535         !! extra MEDUSA-2 tracers 
    452536         DO jj = 2,jpjm1 
    453537            DO ji = 2,jpim1 
     
    456540                  zdtc(ji,jj) = max(0.,trn(ji,jj,jk,jpdtc)) 
    457541                  !! dissolved inorganic carbon 
    458                   zdic(ji,jj) = max(0.,trn(ji,jj,jk,jpdic)) 
     542                  zdic(ji,jj) = trn(ji,jj,jk,jpdic) 
    459543                  !! alkalinity 
    460                   zalk(ji,jj) = max(0.,trn(ji,jj,jk,jpalk)) 
     544                  zalk(ji,jj) = trn(ji,jj,jk,jpalk) 
    461545                  !! oxygen 
    462546                  zoxy(ji,jj) = max(0.,trn(ji,jj,jk,jpoxy)) 
     
    470554                  ztmp(ji,jj) = tsn(ji,jj,jk,jp_tem) 
    471555                  zsal(ji,jj) = tsn(ji,jj,jk,jp_sal) 
    472                   !! 
    473              !! AXY (28/02/14): check input fields 
    474                   if (ztmp(ji,jj) .lt. -3.0 .or. ztmp(ji,jj) .gt. 40.0 ) then 
    475                      IF(lwp) WRITE(numout,*)                                 & 
    476                         ' trc_bio_medusa: T WARNING 2D, ',                   & 
    477                         tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem),          & 
    478                         ' at (', ji, ',', jj, ',', jk, ') at time', kt 
    479            IF(lwp) WRITE(numout,*)                                 & 
    480                         ' trc_bio_medusa: T SWITCHING 2D, ',                 & 
    481                         tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    482                      !! temperatur 
    483                      ztmp(ji,jj) = tsb(ji,jj,jk,jp_tem) 
    484                   endif 
    485                   if (zsal(ji,jj) .lt. 0.0 .or. zsal(ji,jj) .gt. 45.0 ) then 
    486                      IF(lwp) WRITE(numout,*)                                 & 
    487                         ' trc_bio_medusa: S WARNING 2D, ',                   & 
    488                         tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal),          & 
    489                         ' at (', ji, ',', jj, ',', jk, ') at time', kt 
    490                   endif 
    491556               ENDIF 
    492557            ENDDO 
    493558         ENDDO 
    494559# else 
     560         !! diagnostic MEDUSA-1 detrital carbon tracer 
    495561         DO jj = 2,jpjm1 
    496562            DO ji = 2,jpim1 
    497                if (tmask(ji,jj,jk) == 1) then 
     563               IF (tmask(ji,jj,jk) == 1) THEN 
    498564                  !! implicit detrital carbon 
    499565                  zdtc(ji,jj) = zdet(ji,jj) * xthetad 
     
    502568         ENDDO 
    503569# endif 
     570 
     571# if defined key_roam 
     572         !! --------------------------------------------- 
     573         !! JPALM -- 14-12-2017 -- Here, before any exeptionnal crazy value is 
     574         !!              removed, we want to tell to the Master Processor that this  
     575         !!              Exceptionnal value did exist 
     576         !! 
     577         Call trc_bio_check(kt, jk) 
     578 
     579         !!================================================================ 
     580    !! AXY (03/11/17): check input fields 
     581         !! tracer values that exceed thresholds can cause carbonate system 
     582         !! failures when passed to MOCSY; temporary temperature excursions 
     583         !! in recent UKESM0.8 runs trigger such failures but are too short 
     584         !! to have physical consequences; this section checks for such 
     585         !! values and amends them using neighbouring values 
     586         !!  
     587         !! the check on temperature here is also carried out at the end of 
     588         !! each model time step and anomalies are reported in the master 
     589         !! ocean.output file; the error reporting below is strictly local 
     590         !! to the relevant ocean.output_XXXX file so will not be visible 
     591         !! unless all processors are reporting output 
     592         !!================================================================ 
     593         !! 
     594         DO jj = 2,jpjm1 
     595            DO ji = 2,jpim1 
     596               if (tmask(ji,jj,jk) == 1) then 
     597                  !! the thresholds for the four tracers are ... 
     598                  IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT. 40.0  ) .OR.   & 
     599                       (zsal(ji,jj) .LE.  0.0) .OR. (zsal(ji,jj) .GT. 50.0  ) .OR.   & 
     600                       (zdic(ji,jj) .LE.  0.0) .OR. (zdic(ji,jj) .GT. 4.0E3 ) .OR.   & 
     601                       (zalk(ji,jj) .LE.  0.0) .OR. (zalk(ji,jj) .GT. 4.0E3 ) ) THEN 
     602                     !! 
     603                     !! all tracer values are reported in the event of any excursion 
     604                     IF (lwp) THEN 
     605                         WRITE(charout,*)  ' Tmp = ', ztmp(ji,jj) 
     606                         WRITE(charout2,*) ' Sal = ', zsal(ji,jj) 
     607                         WRITE(charout3,*) ' DIC = ', zdic(ji,jj) 
     608                         WRITE(charout4,*) ' Alk = ', zalk(ji,jj) 
     609                         WRITE(charout5,*) mig(ji), mjg(jj), jk, kt  
     610                         CALL ctl_warn( 'trc_bio_medusa: carbonate chemistry WARNING:',  & 
     611                            TRIM(charout),TRIM(charout2),TRIM(charout3),TRIM(charout4),  &  
     612                            'at i, j, k, kt:', TRIM(charout5) ) 
     613                     ENDIF 
     614                     !! 
     615                     !! Detect, report and correct tracer excursions 
     616                     IF ( (ztmp(ji,jj) .LT. -3.0) .OR. (ztmp(ji,jj) .GT.  40.0) )             & 
     617                        CALL trc_bio_exceptionnal_fix(                                         & 
     618                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_tem), tmask(ji-1:ji+1,jj-1:jj+1,jk),    & 
     619                        'Tmp', -3.0, 40.0, ztmp(ji,jj) ) 
     620                     !! 
     621                     IF ( (zsal(ji,jj) .LE. 0.0) .OR. (zsal(ji,jj) .GT.  50.0)  )             & 
     622                        CALL trc_bio_exceptionnal_fix(                                         & 
     623                        tsn(ji-1:ji+1,jj-1:jj+1,jk,jp_sal), tmask(ji-1:ji+1,jj-1:jj+1,jk),    & 
     624                        'Sal', 1.0, 50.0, zsal(ji,jj) ) 
     625                     !! 
     626                     IF ( (zdic(ji,jj) .LE. 0.0) .OR. (zdic(ji,jj) .GT. 4.0E3)  )             & 
     627                        CALL trc_bio_exceptionnal_fix(                                         & 
     628                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpdic), tmask(ji-1:ji+1,jj-1:jj+1,jk),     & 
     629                        'DIC', 100.0, 4.0E3, zdic(ji,jj) ) 
     630                     !! 
     631                     IF ( (zalk(ji,jj) .LE. 0.0) .OR. (zalk(ji,jj) .GT. 4.0E3)  )             & 
     632                        CALL trc_bio_exceptionnal_fix(                                         & 
     633                        trn(ji-1:ji+1,jj-1:jj+1,jk,jpalk), tmask(ji-1:ji+1,jj-1:jj+1,jk),     & 
     634                        'Alk', 100.0, 4.0E3, zalk(ji,jj) ) 
     635                  ENDIF 
     636               ENDIF 
     637            ENDDO 
     638         ENDDO 
     639# endif 
     640 
    504641# if defined key_debug_medusa 
    505642         DO jj = 2,jpjm1 
     
    657794   END SUBROUTINE trc_bio_medusa 
    658795 
     796 
     797 
     798   SUBROUTINE trc_bio_exceptionnal_fix(tiny_var, tiny_mask, var_nm, mini, maxi, varout) 
     799      !! JPALM (27/10/17): This function is called only when abnormal values that  
     800      !! could break the model's carbonate system are fed to MEDUSA 
     801      !!  
     802      !! The function calculates an average tracer value based on the 3x3 cell 
     803      !! neighbourhood around the abnormal cell, and reports this back 
     804      !!  
     805      !! Tracer array values are not modified, but MEDUSA uses "corrected" values 
     806      !! in its calculations 
     807      !! 
     808      !! temporary variables 
     809      REAL(wp), INTENT( in ), DIMENSION(3,3) :: tiny_var, tiny_mask 
     810      CHARACTER(3), INTENT( in )            :: var_nm 
     811      REAL(wp), INTENT( in )                 :: mini, maxi 
     812      REAL(wp), INTENT( out )                :: varout 
     813      REAL(wp)      :: sumtsn, tsnavg 
     814      INTEGER       :: summask 
     815      CHARACTER(25) :: charout1, charout2 
     816      CHARACTER(60) :: charout3, charout4 
     817      INTEGER       :: ii,ij 
     818     
     819      !! point to the center of the 3*3 zoom-grid, to check around 
     820      ii = 2 
     821      ij = 2 
     822      !! Print surounding values to check if isolated Crazy value or  
     823      !! General error 
     824      IF(lwp) THEN  
     825          WRITE(numout,*)                                 & 
     826            '----------------------------------------------------------------------' 
     827          WRITE(numout,*)                                 & 
     828            'trc_bio_medusa: 3x3 neighbourhood surrounding abnormal ', TRIM(var_nm) 
     829          WRITE(numout,9100)                              & 
     830            3, tiny_var(ii-1,ij+1), tiny_var(ii  ,ij+1), tiny_var(ii+1,ij+1) 
     831          WRITE(numout,9100)                              & 
     832            2, tiny_var(ii-1,ij  ), tiny_var(ii  ,ij  ), tiny_var(ii+1,ij  ) 
     833          WRITE(numout,9100)                              & 
     834            1, tiny_var(ii-1,ij-1), tiny_var(ii  ,ij-1), tiny_var(ii+1,ij-1) 
     835          WRITE(numout,*)                                 & 
     836            'trc_bio_medusa: 3x3 land-sea neighbourhood, tmask' 
     837          WRITE(numout,9100)                              & 
     838            3, tiny_mask(ii-1,ij+1), tiny_mask(ii  ,ij+1), tiny_mask(ii+1,ij+1) 
     839          WRITE(numout,9100)                              & 
     840            2, tiny_mask(ii-1,ij  ), tiny_mask(ii  ,ij  ), tiny_mask(ii+1,ij  ) 
     841          WRITE(numout,9100)                              & 
     842            1, tiny_mask(ii-1,ij-1), tiny_mask(ii  ,ij-1), tiny_mask(ii+1,ij-1) 
     843      ENDIF 
     844      !! Correct out of range values 
     845      sumtsn = ( tiny_mask(ii-1,ij+1) * tiny_var(ii-1,ij+1) ) +  & 
     846               ( tiny_mask(ii  ,ij+1) * tiny_var(ii  ,ij+1) ) +  & 
     847               ( tiny_mask(ii+1,ij+1) * tiny_var(ii+1,ij+1) ) +  & 
     848               ( tiny_mask(ii-1,ij  ) * tiny_var(ii-1,ij  ) ) +  & 
     849               ( tiny_mask(ii+1,ij  ) * tiny_var(ii+1,ij  ) ) +  & 
     850               ( tiny_mask(ii-1,ij-1) * tiny_var(ii-1,ij-1) ) +  & 
     851               ( tiny_mask(ii  ,ij-1) * tiny_var(ii  ,ij-1) ) +  & 
     852               ( tiny_mask(ii+1,ij-1) * tiny_var(ii+1,ij-1) ) 
     853      !! 
     854      summask = tiny_mask(ii-1,ij+1) + tiny_mask(ii  ,ij+1)   +  & 
     855                tiny_mask(ii+1,ij+1) + tiny_mask(ii-1,ij  )   +  & 
     856                tiny_mask(ii+1,ij  ) + tiny_mask(ii-1,ij-1)   +  & 
     857                tiny_mask(ii  ,ij-1) + tiny_mask(ii+1,ij-1) 
     858      !! 
     859      IF ( summask .GT. 0 ) THEN 
     860         tsnavg = ( sumtsn / summask ) 
     861         varout = MAX( MIN( maxi, tsnavg), mini ) 
     862      ELSE    
     863         IF (ztmp(ii,ij) .LT. mini )  varout = mini 
     864         IF (ztmp(ii,ij) .GT. maxi )  varout = maxi 
     865      ENDIF 
     866      !! 
     867      IF (lwp) THEN  
     868          WRITE(charout1,9200) tiny_var(ii,ij) 
     869          WRITE(charout2,9200) varout 
     870          WRITE(charout3,*) ' ', charout1, ' -> ', charout2 
     871          WRITE(charout4,*) ' Tracer: ', trim(var_nm) 
     872      !! 
     873          WRITE(numout,*) 'trc_bio_medusa: ** EXCEPTIONAL VALUE SWITCHING **' 
     874          WRITE(numout,*) charout4  
     875          WRITE(numout,*) charout3 
     876          WRITE(numout,*) '----------------------------------------------------------------------' 
     877          WRITE(numout,*) ' ' 
     878      ENDIF 
     879 
     8809100  FORMAT('Row:', i1, '  ', e16.6, ' ', e16.6, ' ', e16.6) 
     8819200  FORMAT(e16.6) 
     882 
     883   END SUBROUTINE trc_bio_exceptionnal_fix  
     884 
     885   SUBROUTINE trc_bio_check(kt, jk) 
     886      !!----------------------------------- 
     887      !! JPALM -- 14-12-2017 -- Still dealing with this micro-boil/carb failure 
     888      !!                     problem. The model is now able to correct a local 
     889      !!                     crazy value. but does it silently. 
     890      !!                     We need to spread the word to the master processor. we 
     891      !!                     don't want the model  to correct values without telling us 
     892      !!                     This module will tell at least when crazy DIC or 
     893      !!                     ALK appears. 
     894      !!------------------------------------- 
     895      REAL(wp)              :: zmax, zmin    ! temporary scalars 
     896      INTEGER               :: ji,jj         ! dummy loop 
     897      INTEGER               :: ii,ij         ! temporary scalars  
     898      INTEGER, DIMENSION(2) :: ilocs         !  
     899      INTEGER, INTENT( in ) :: kt, jk 
     900      !! 
     901      !!========================== 
     902      !! DIC Check 
     903      zmax =  -5.0  ! arbitrary  low maximum value 
     904      zmin =  4.0E4  ! arbitrary high minimum value 
     905      DO jj = 2, jpjm1 
     906         DO ji = 2,jpim1 
     907            IF( tmask(ji,jj,1) == 1) THEN 
     908               zmax = MAX(zmax,zdic(ji,jj))     ! find local maximum 
     909               zmin = MIN(zmin,zdic(ji,jj))     ! find local minimum 
     910            ENDIF 
     911         END DO 
     912      END DO 
     913      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain 
     914      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain 
     915      ! 
     916      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem 
     917         IF (lk_mpp) THEN 
     918            CALL mpp_maxloc ( zdic(:,:),tmask(:,:,1), zmax, ii,ij ) 
     919         ELSE 
     920            ilocs = MAXLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 
     921            ii = ilocs(1) + nimpp - 1 
     922            ij = ilocs(2) + njmpp - 1 
     923         ENDIF 
     924         ! 
     925         IF(lwp) THEN 
     926            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     927            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration > 4000 ' 
     928            WRITE(numout,9600) kt, zmax, ii, ij, jk 
     929            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     930         ENDIF 
     931      ENDIF 
     932      ! 
     933      IF( zmin .LE. 0.0) THEN  ! we've got a problem 
     934         IF (lk_mpp) THEN 
     935            CALL mpp_minloc ( zdic(:,:),tmask(:,:,1), zmin, ii,ij ) 
     936         ELSE 
     937            ilocs = MINLOC( zdic(:,:), mask = tmask(:,:,1) == 1. ) 
     938            ii = ilocs(1) + nimpp - 1 
     939            ij = ilocs(2) + njmpp - 1 
     940         ENDIF 
     941         ! 
     942         IF(lwp) THEN 
     943            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     944            WRITE(numout,*) 'trc_bio:tracer anomaly: DIC concentration <= 0 ' 
     945            WRITE(numout,9700) kt, zmin, ii, ij, jk 
     946            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     947         ENDIF 
     948      ENDIF 
     949      !! 
     950      !!========================== 
     951      !! ALKALINITY Check 
     952      zmax =  -5.0  ! arbitrary  low maximum value 
     953      zmin =  4.0E4  ! arbitrary high minimum value 
     954      DO jj = 2, jpjm1 
     955         DO ji = 2,jpim1 
     956            IF( tmask(ji,jj,1) == 1) THEN 
     957               zmax = MAX(zmax,zalk(ji,jj))     ! find local maximum 
     958               zmin = MIN(zmin,zalk(ji,jj))     ! find local minimum 
     959            ENDIF 
     960         END DO 
     961      END DO 
     962      IF( lk_mpp )   CALL mpp_max( zmax )       ! max over the global domain 
     963      IF( lk_mpp )   CALL mpp_min( zmin )       ! min over the global domain 
     964      ! 
     965      IF( zmax .GT. 4.0E3) THEN  ! we've got a problem 
     966         IF (lk_mpp) THEN 
     967            CALL mpp_maxloc ( zalk(:,:),tmask(:,:,1), zmax, ii,ij ) 
     968         ELSE 
     969            ilocs = MAXLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 
     970            ii = ilocs(1) + nimpp - 1 
     971            ij = ilocs(2) + njmpp - 1 
     972         ENDIF 
     973         ! 
     974         IF(lwp) THEN 
     975            WRITE(numout,*) 'trc_bio:tracer anomaly: *****     WARNING    *****' 
     976            WRITE(numout,*) 'trc_bio:tracer anomaly: ALK concentration > 4000 ' 
     977            WRITE(numout,9800) kt, zmax, ii, ij, jk 
     978            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     979         ENDIF 
     980      ENDIF 
     981      ! 
     982      IF( zmin .LE. 0.0) THEN  ! we've got a problem 
     983         IF (lk_mpp) THEN 
     984            CALL mpp_minloc ( zalk(:,:),tmask(:,:,1), zmin, ii,ij ) 
     985         ELSE 
     986            ilocs = MINLOC( zalk(:,:), mask = tmask(:,:,1) == 1. ) 
     987            ii = ilocs(1) + nimpp - 1 
     988            ij = ilocs(2) + njmpp - 1 
     989         ENDIF 
     990         ! 
     991         IF(lwp) THEN 
     992            WRITE(numout,*) 'trc_bio:tracer anomaly: *****    WARNING     *****' 
     993            WRITE(numout,*) 'trc_bio:tracer anomaly:  ALK concentration <= 0 ' 
     994            WRITE(numout,9900) kt, zmin, ii, ij, jk 
     995            WRITE(numout,*) 'trc_bio:tracer anomaly: ***** END OF WARNING *****' 
     996         ENDIF 
     997      ENDIF 
     998 
     999 
     10009600  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max DIC: ',f16.10,', i j k: ',3i5) 
     10019700  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min DIC: ',f16.10,', i j k: ',3i5) 
     10029800  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' max ALK: ',f16.10,', i j k: ',3i5) 
     10039900  FORMAT ('trc_bio:tracer anomaly: kt=',i6,' min ALK: ',f16.10,', i j k: ',3i5) 
     1004 
     1005   END SUBROUTINE trc_bio_check 
     1006 
     1007 
    6591008#else 
    6601009   !!===================================================================== 
     
    6631012CONTAINS 
    6641013   SUBROUTINE trc_bio_medusa( kt )                   ! Empty routine 
     1014      IMPLICIT NONE 
    6651015      INTEGER, INTENT( in ) ::   kt 
    6661016      WRITE(*,*) 'trc_bio_medusa: You should not have seen this print! error?', kt 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcdms_medusa.F90

    r8132 r10232  
    113113! 
    114114! AXY (13/03/15): Anderson et al. (2001) 
     115!! JPALM --19-12-2017-- Tunable through the namelist 
     116!!                      within dmsmin - dmscut - dmsslp 
    115117        Jterm = xqsr + 1.0e-6 
    116118        !! this next line makes a hard-coded assumption about the  
     
    120122        Qterm = xdin / (xdin + 0.5) 
    121123        fq1 = log10(CHL * Jterm * Qterm) 
    122         if (fq1 > 1.72) then 
    123            dms_andr = (8.24 * (fq1 - 1.72)) + 2.29 
    124         else 
    125            dms_andr = 2.29 
     124        if (fq1 > dmscut) then 
     125           dms_andr = (dmsslp * (fq1 - dmscut)) + dmsmin 
     126        else 
     127           dms_andr = dmsmin 
    126128        endif 
    127129! 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90

    r8147 r10232  
    271271      zn_dms_srf(:,:)  = 0.0 
    272272      za_dms_srf(:,:)  = 0.0 
    273       zn_chl_srf(:,:)  = 2.0E-8 !! Chl srf 
     273      zn_chl_srf(:,:)  = 2.0E-8 !! Chl cpl - set first as surf 
    274274      !! 
    275275      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: DMS and CO2 flux (UKESM) initialised to zero' 
     
    278278         CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable 
    279279         DMS_out_cpl(:,:)     =  zn_dms_srf(:,:)   !! Coupling variable 
    280          chloro_out_cpl(:,:)  =  zn_chl_srf(:,:)   !! Coupling variable 
     280         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:) * scl_chl   !! Coupling variable 
    281281      END IF 
    282282      !! 
     
    324324      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    325325      IF(lwp) CALL flush(numout) 
     326      !! 
     327      !!---------------------------------------------------------------------- 
     328      !! JPALM (23-01-2017): new way to initialize CO2-atm for cmip6  
     329      !!                     initially done in trcsms_medusa 
     330      !!---------------------------------------------------------------------- 
     331      !!  
     332      IF( ( .NOT.lk_oasis ) .AND. ( .NOT.lk_pi_co2 ) .AND. ( xobs_xco2a <= 0.0 ) ) THEN 
     333         IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisating atm CO2 record' 
     334         CALL trc_ini_medusa_co2atm 
     335      ENDIF 
    326336 
    327337   END SUBROUTINE trc_ini_medusa 
     
    480490   END SUBROUTINE trc_ini_medusa_river 
    481491    
     492   SUBROUTINE trc_ini_medusa_co2atm 
     493      !!---------------------------------------------------------------------- 
     494      !!                     ***  trc_ini_medusa_co2atm  ***   
     495      !! 
     496      !! ** Purpose :   initialization atmospheric co2 record 
     497      !! 
     498      !! ** Method  : - Read the xco2 file 
     499      !!---------------------------------------------------------------------- 
     500      INTEGER                       ::  jn, jm, io, ierr, inum, iostatus 
     501      INTEGER, PARAMETER            ::  iskip = 4   ! number of 1st descriptor lines 
     502      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   zyy !: xCO2 record years 
     503      CHARACTER (len=10), PARAMETER ::  clname = 'xco2.atm'  !! atm CO2 record file 
     504      !!---------------------------------------------------------------------- 
     505 
     506      IF(lwp) WRITE(numout,*) 
     507      IF(lwp) WRITE(numout,*) ' trc_ini_medusa_co2atm: initialisation of atm CO2 historical record' 
     508      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 
     509 
     510 
     511      IF(lwp) WRITE(numout,*) 'read of formatted file xco2.atm' 
     512 
     513      CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
     514      !!! 
     515      ! -Compute the number of year in the file 
     516      ! -File starts in co2_yinit, jn represents the record number in the file. 
     517      ! -Remove the file head (iskip lines) to jn 
     518      ! -The year is jn + yinit - 1  
     519      !! Determine the number of lines in xCO2 input file 
     520      iostatus = 0 
     521      jn = 1 
     522      DO WHILE ( iostatus == 0 ) 
     523        READ(inum,'(1x)', IOSTAT=iostatus, END=100) 
     524        jn = jn + 1 
     525      ENDDO 
     526      IF( iostatus .NE. 0 ) THEN 
     527        !! Error while reading xCO2 input file  
     528        CALL ctl_stop('trc_ini_medusa_co2atm: & 
     529                      & Error on the 1st reading of xco2.atm') 
     530        RETURN 
     531      ENDIF 
     532 100  co2_rec = jn - 1 - iskip 
     533      IF ( lwp) WRITE(numout,*) '    ', co2_rec ,' years read in the file' 
     534      !                                ! Allocate CO2 hist arrays 
     535      ierr = 0  
     536      ALLOCATE( hist_pco2(co2_rec),zyy(co2_rec), STAT=ierr ) 
     537      IF( ierr > 0 ) THEN 
     538         CALL ctl_stop( 'trc_ini_medusa_co2atm: unable to allocate  array' )   
     539         RETURN 
     540      ENDIF 
     541 
     542      REWIND(inum) 
     543 
     544      DO jm = 1, iskip        ! Skip over 1st six descriptor lines 
     545         READ(inum,'(1x)') 
     546      END DO 
     547      ! file starts in 1931 do jn represent the year in the century.jhh 
     548      ! Read file till the end 
     549      ! allocate start and end year of the file 
     550      DO jn = 1, co2_rec 
     551        READ(inum,'(F6.1,F12.7)', IOSTAT=io) zyy(jn), hist_pco2(jn) 
     552        IF( io .NE. 0 ) THEN 
     553          !! Error while reading xCO2 input file  
     554          CALL ctl_stop('trc_ini_medusa_co2atm: & 
     555                        & Error on the 2nd reading of xco2.atm') 
     556          RETURN 
     557        ENDIF 
     558 
     559        IF(jn==1) co2_yinit = zyy(jn) 
     560      END DO 
     561      co2_yend = co2_yinit + real(co2_rec) - 1. 
     562 
     563      IF(lwp) THEN        ! Control print 
     564         WRITE(numout,*) 
     565         WRITE(numout,*) 'CO2 hist start year: ', co2_yinit 
     566         WRITE(numout,*) 'CO2 hist end   year: ', co2_yend 
     567         WRITE(numout,*) ' Year   xCO2 atm ' 
     568         DO jn = 1, co2_rec 
     569            WRITE(numout, '(F6.1,F12.7)') zyy(jn), hist_pco2(jn) 
     570         END DO 
     571      ENDIF 
     572 
     573   END SUBROUTINE trc_ini_medusa_co2atm 
     574 
     575 
    482576#else 
    483577   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90

    r8442 r10232  
    2424   USE sms_medusa      ! sms trends 
    2525   USE iom             ! I/O manager 
     26   USE sbc_oce, ONLY: lk_oasis 
    2627   !!USE trc_nam_dia     ! JPALM 13-11-2015 -- if iom_use for diag 
    2728 
     
    8586      &  jriver_n,jriver_si,jriver_c,jriver_alk,jriver_dep,   & 
    8687      &  xsdiss,                                              & 
    87       &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model 
     88      &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model,   & 
     89      &  scl_chl, chl_out, dmsmin, dmscut, dmsslp,            & 
     90      &  ln_foam_medusa 
    8891#if defined key_roam 
    8992      NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        & 
     93      &    xobs_xco2a,                                        & 
    9094      &    xthetarem,xo2min  
    9195#endif 
     
    125129      ! 1.4 namelist natbio : biological parameters 
    126130      ! ------------------------------------------- 
     131      !! Note: the default values below will all be overwritten by the 
     132      !!       input in the namelist natbio. 
    127133       
    128134      xxi         = 0. 
     
    246252      jdms_input  = 0 
    247253      jdms_model  = 0 
     254      scl_chl     = 1. 
     255      chl_out     = 1 
     256      dmsmin      = 2.29 !! Anderson DMS default 
     257      dmscut      = 1.72 !! Anderson DMS default 
     258      dmsslp      = 8.24 !! Anderson DMS default 
     259!! 
     260      ln_foam_medusa = .FALSE. 
    248261             
    249262      !REWIND(numnatm) 
     
    399412!!       jdms_model  :  choice of DMS model passed to atmosphere 
    400413!!                      1 = ANDR, 2 = SIMO, 3 = ARAN, 4 = HALL, 5 = ANDM 
    401 !! 
     414!!       dmsmin      : DMS minimum value for DMS Anderson (ANDR) sheme ONLY 
     415!!       dmscut      : DMS cutoff value for DMS Anderson (ANDR) sheme ONLY 
     416!!       dmsslp      : DMS slope value for DMS Anderson (ANDR) sheme ONLY 
     417!! UKESM1 - exported Chl to UM 
     418!!       scl_chl     : scaling factor to tune the chl field sent to the UM 
     419!!       chl_out     : select the chl field to send at the UM: 
     420!!                     1- Surf Chl ; 2- MLD Chl  
     421!! 
     422!! FOAM - observation operator and data assimilation 
     423!!       ln_foam_medusa : calculate required diagnostics 
     424 
    402425      IF(lwp) THEN 
    403426!! 
     
    909932!! 
    910933!! UKESM1 - new diagnostics  !! Jpalm; AXY (08/07/15) 
    911          WRITE(numout,*) '=== UKESM1-related parameters' 
    912          WRITE(numout,*)     & 
    913          &   ' include DMS diagnostic?,                                   jdms        = ', jdms 
    914          if (jdms_input .eq. 0) then 
    915             WRITE(numout,*)     & 
    916             &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous' 
    917          else 
    918             WRITE(numout,*)     & 
    919             &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average' 
    920          endif 
    921     if (jdms_model .eq. 1) then 
    922             WRITE(numout,*)     & 
    923             &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)' 
    924     elseif (jdms_model .eq. 2) then 
    925             WRITE(numout,*)     & 
    926             &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)' 
    927     elseif (jdms_model .eq. 3) then 
    928             WRITE(numout,*)     & 
    929             &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)' 
    930     elseif (jdms_model .eq. 4) then 
    931             WRITE(numout,*)     & 
    932             &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)' 
    933     elseif (jdms_model .eq. 5) then 
    934             WRITE(numout,*)     & 
    935             &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)' 
    936          endif 
     934         WRITE(numout,*) '=== UKESM1-related parameters ===' 
     935         WRITE(numout,*) ' ---- --- ---' 
     936 
     937         IF (lk_oasis) THEN 
     938            WRITE(numout,*) '=== UKESM1 --  coupled DMS to the atmosphere' 
     939            WRITE(numout,*)     & 
     940            &   ' include DMS diagnostic?,                                   jdms        = ', jdms 
     941            if (jdms_input .eq. 0) then 
     942               WRITE(numout,*)     & 
     943               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = instantaneous' 
     944            else 
     945               WRITE(numout,*)     & 
     946               &   ' use instant (0) or diel-avg (1) inputs,                    jdms_input  = diel-average' 
     947            endif 
     948          if (jdms_model .eq. 1) then 
     949               WRITE(numout,*)     & 
     950               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001)' 
     951       elseif (jdms_model .eq. 2) then 
     952               WRITE(numout,*)     & 
     953               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Simo & Dachs (2002)' 
     954       elseif (jdms_model .eq. 3) then 
     955               WRITE(numout,*)     & 
     956               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Aranami & Tsunogai (2004)' 
     957       elseif (jdms_model .eq. 4) then 
     958               WRITE(numout,*)     & 
     959               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Halloran et al. (2010)' 
     960       elseif (jdms_model .eq. 5) then 
     961               WRITE(numout,*)     & 
     962               &   ' choice of DMS model passed to atmosphere,                  jdms_model  = Anderson et al. (2001; modified)' 
     963            endif 
     964            if (jdms_model .eq. 1) then 
     965               WRITE(numout,*)     & 
     966               &   ' Anderson DMS model tuned parameters:                       DMS minimum = ',dmsmin,'. -- Default = 2.29 ' 
     967               WRITE(numout,*)     & 
     968               &   ' Anderson DMS model tuned parameters:                       DMS cutoff  = ',dmscut,'. -- Default = 1.72 ' 
     969               WRITE(numout,*)     & 
     970               &   ' Anderson DMS model tuned parameters:                       DMS slope   = ',dmsslp,'. -- Default = 8.24 ' 
     971            endif 
     972 
     973            WRITE(numout,*) '=== UKESM1 --  coupled Chl to the atmosphere' 
     974            WRITE(numout,*)        & 
     975               &   ' Scaling factor to export tuned Chl to the atmosphere       scl_chl  = ', scl_chl 
     976            IF (chl_out .eq. 1) THEN 
     977               WRITE(numout,*)        & 
     978               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = Surface Chl field ' 
     979            ELSEIF (chl_out .eq. 2) THEN 
     980               WRITE(numout,*)        & 
     981               &   ' Chl field to be scaled and sent to the atmosphere:         chl_out  = MLD Chl field ' 
     982            ENDIF 
     983         ENDIF ! IF lk_oasis=true 
     984!! FOAM 
     985         WRITE(numout,*) '=== FOAM-related parameters' 
     986         WRITE(numout,*)     & 
     987         &   ' calculate diagnostics for data assimilation,            ln_foam_medusa = ', ln_foam_medusa 
    937988!! 
    938989      ENDIF 
     
    9901041      xthetarem = 0. 
    9911042      xo2min    = 0. 
     1043      xobs_xco2a = 0. 
    9921044 
    9931045      !READ(numnatm,natroam) 
     
    10091061!!       xthetarem :  oxygen consumption by carbon remineralisation 
    10101062!!       xo2min    :  oxygen minimum concentration 
     1063!!       xobs_xco2a : observed atmospheric xCO2 (not used if <= 0) 
    10111064 
    10121065      IF(lwp) THEN 
     
    10261079          WRITE(numout,*)     & 
    10271080          &   ' oxygen minimum concentration                               xo2min      = ', xo2min 
     1081          WRITE(numout,*)     & 
     1082          &   ' observed atmospheric xCO2 (not used if <= 0)               xobs_xco2a  = ', xobs_xco2a 
    10281083       ENDIF 
    10291084 
     
    20532108          med_diag%OCN_DPCO2%dgsave = .FALSE. 
    20542109      ENDIF 
    2055       !! 
     2110      !! UKESM additional 
     2111      IF  (iom_use("CHL_MLD")) THEN  
     2112          med_diag%CHL_MLD%dgsave = .TRUE. 
     2113      ELSE  
     2114          med_diag%CHL_MLD%dgsave = .FALSE. 
     2115      ENDIF 
     2116      IF  (iom_use("CHL_CPL")) THEN  
     2117          med_diag%CHL_CPL%dgsave = .TRUE. 
     2118      ELSE  
     2119          med_diag%CHL_CPL%dgsave = .FALSE. 
     2120      ENDIF 
     2121      !! 3D 
    20562122      IF  (iom_use("TPP3")) THEN  
    20572123          med_diag%TPP3%dgsave = .TRUE. 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcsms_medusa.F90

    r8074 r10232  
    88   !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA 
    99   !!              -   !  2010-03  (A. Yool) updated for branch inclusion 
     10   !!              -   !  2017-08  (A. Yool) amend for slow detritus bug 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_medusa 
     
    2223   USE trcsed_medusa 
    2324   USE trcavg_medusa 
     25   !! for SMS trends 
     26   USE par_medusa,    ONLY: jp_msa0, jp_msa1, jp_medusa 
     27   USE par_oce,       ONLY: jpi, jpj, jpk 
     28   USE trd_oce,       ONLY: jptra_sms, l_trdtrc 
     29   USE trdtrc 
    2430 
    2531 
     
    4652      !!---------------------------------------------------------------------- 
    4753      INTEGER, INTENT(in) :: kt   ! ocean time-step index 
     54      !! Loop variables 
     55      INTEGER :: jn 
     56      !! trend temporary array: 
     57      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrmed 
     58 
    4859 
    4960# if defined key_debug_medusa 
     
    5768       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    5869      ENDIF 
     70 
     71      !! MEDUSA SMS trends: 
     72      IF( l_trdtrc ) THEN 
     73          CALL wrk_alloc( jpi, jpj, jpk, jp_medusa, ztrmed ) 
     74          ztrmed(:,:,:,:)=0.0   
     75          DO jn = 1, jp_medusa 
     76            ztrmed(:,:,:,jn) = tra(:,:,:,jp_msa0 + jn - 1) 
     77          END DO 
     78      END IF 
    5979 
    6080      CALL trc_avg_medusa( kt ) ! rolling average module 
     
    88108#  endif 
    89109       
    90       CALL trc_sed_medusa( kt ) ! sedimentation model 
    91 #  if defined key_debug_medusa 
    92       IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa' 
    93       CALL flush(numout) 
    94 #  endif 
     110!! AXY (08/08/2017): remove call to buggy subroutine (now handled by detritus.F90) 
     111!!       CALL trc_sed_medusa( kt ) ! sedimentation model 
     112!! #  if defined key_debug_medusa 
     113!!       IF(lwp) WRITE(numout,*) ' MEDUSA done trc_sed_medusa' 
     114!!       CALL flush(numout) 
     115!! #  endif 
    95116# endif 
     117 
     118      !! MEDUSA SMS trends: 
     119      IF( l_trdtrc ) THEN 
     120          DO jn = 1, jp_medusa 
     121            ztrmed(:,:,:,jn) = tra(:,:,:,jp_msa0 + jn - 1)-ztrmed(:,:,:,jn) 
     122            CALL trd_trc( ztrmed(:,:,:,jn), jn, jptra_sms, kt )   ! save trends 
     123          END DO 
     124          CALL wrk_dealloc( jpi, jpj, jpk, jp_medusa, ztrmed ) 
     125      END IF 
     126 
    96127 
    97128   END SUBROUTINE trc_sms_medusa 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_bgcupdates/NEMOGCM/NEMO/TOP_SRC/MEDUSA/zooplankton.F90

    r8441 r10232  
    4141                                   idf, idfval,                          & 
    4242                                   zdet, zdtc, zphd, zphn, zzme, zzmi 
    43       USE dom_oce,           ONLY: e3t_0, e3t_n, tmask 
     43      USE dom_oce,           ONLY: e3t_0, tmask 
     44#if defined key_vvl 
     45      USE dom_oce,           ONLY: e3t_n 
     46#endif 
    4447      USE par_kind,          ONLY: wp 
    4548      USE in_out_manager,    ONLY: lwp, numout 
Note: See TracChangeset for help on using the changeset viewer.