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 8657 – NEMO

Changeset 8657


Ignore:
Timestamp:
2017-10-25T14:46:18+02:00 (7 years ago)
Author:
jpalmier
Message:

update the branch to match last GO6 changes

Location:
branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM
Files:
2 deleted
39 edited
18 copied

Legend:

Unmodified
Added
Removed

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/CONFIG/SHARED/field_def_bgc.xml

    r8308 r8657  
    448448       <field id= "DMS_HALL"   long_name="DMS Surface Concentration, Halloran"       unit="nmol/L"      /> 
    449449       <field id= "DMS_ANDM"   long_name="DMS Surface Concentration, Anderson modif" unit="nmol/L"      /> 
     450       <field id= "CHL_MLD"    long_name="MLD averaged Chlorophyll"                  unit="mg Chl/m3"   /> 
    450451       <field id= "ATM_XCO2"   long_name="Atmospheric xCO2"                          unit="ppm"         /> 
    451452       <field id= "OCN_FCO2"   long_name="Surface ocean fCO2"                        unit="uatm"        /> 
    … …  
    784785      <field field_ref= "CO2STARAIR" name="CO2STARAIR" /> 
    785786      <field field_ref= "OCN_DPCO2"  name="OCN_DPCO2"  /> 
     787      <field field_ref= "CHL_MLD"    name="CHL_MLD"    /> 
    786788    </field_group> 
    787789 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8280 r8657  
    971971                           !        = 0  constant 10 m length scale 
    972972                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
     973   rn_c        =   0.8     !  Default value only used when nn_htau = 2 (typically never!) 
    973974/ 
    974975!------------------------------------------------------------------------ 
    … …  
    12251226   ln_s3d     = .false.    ! Logical switch for S profile observations 
    12261227   ln_ena     = .false.    ! Logical switch for ENACT insitu data set 
    1227    ln_cor     = .false.    ! Logical switch for Coriolis insitu data set 
     1228   !                       !     ln_cor                  Logical switch for Coriolis insitu data set 
    12281229   ln_profb   = .false.    ! Logical switch for feedback insitu data set 
    12291230   ln_sla     = .false.    ! Logical switch for SLA observations 
     1231 
    12301232   ln_sladt   = .false.    ! Logical switch for AVISO SLA data 
     1233 
    12311234   ln_slafb   = .false.    ! Logical switch for feedback SLA data 
    1232    ln_ssh     = .false.    ! Logical switch for SSH observations 
    1233    ln_sst     = .false.    ! Logical switch for SST observations 
    1234    ln_reysst  = .false.    ! Logical switch for Reynolds observations 
    1235    ln_ghrsst  = .false.    ! Logical switch for GHRSST observations 
     1235                           !     ln_ssh                  Logical switch for SSH observations 
     1236 
     1237   ln_sst     = .false.     ! Logical switch for SST observations 
     1238   ln_reysst  = .false.     !     ln_reysst               Logical switch for Reynolds observations 
     1239   ln_ghrsst  = .false.    !     ln_ghrsst               Logical switch for GHRSST observations 
     1240 
    12361241   ln_sstfb   = .false.    ! Logical switch for feedback SST data 
    1237    ln_sss     = .false.    ! Logical switch for SSS observations 
     1242                           !     ln_sss                  Logical switch for SSS observations 
    12381243   ln_seaice  = .false.    ! Logical switch for Sea Ice observations 
    1239    ln_vel3d   = .false.    ! Logical switch for velocity observations 
    1240    ln_velavcur= .false     ! Logical switch for velocity daily av. cur. 
    1241    ln_velhrcur= .false     ! Logical switch for velocity high freq. cur. 
    1242    ln_velavadcp = .false.  ! Logical switch for velocity daily av. ADCP 
    1243    ln_velhradcp = .false.  ! Logical switch for velocity high freq. ADCP 
    1244    ln_velfb   = .false.    ! Logical switch for feedback velocity data 
    1245    ln_grid_global = .false. ! Global distribtion of observations 
    1246    ln_grid_search_lookup = .false. !  Logical switch for obs grid search w/lookup table 
    1247    grid_search_file = 'grid_search'  !  Grid search lookup file header 
    1248 ! All of the *files* variables below are arrays. Use namelist_cfg to add more files 
    1249    enactfiles = 'enact.nc' !  ENACT input observation file names (specify full array in namelist_cfg) 
    1250    coriofiles = 'corio.nc' !  Coriolis input observation file name 
    1251    profbfiles = 'profiles_01.nc' ! Profile feedback input observation file name 
    1252    ln_profb_enatim = .false !        Enact feedback input time setting switch 
    1253    slafilesact = 'sla_act.nc' !  Active SLA input observation file names 
    1254    slafilespas = 'sla_pass.nc' ! Passive SLA input observation file names 
    1255    slafbfiles = 'sla_01.nc' ! slafbfiles: Feedback SLA input observation file names 
    1256    sstfiles = 'ghrsst.nc'   ! GHRSST input observation file names 
    1257    sstfbfiles = 'sst_01.nc' ! Feedback SST input observation file names 
    1258    seaicefiles = 'seaice_01.nc' ! Sea Ice input observation file names 
    1259    velavcurfiles = 'velavcurfile.nc'  ! Vel. cur. daily av. input file name 
    1260    velhrcurfiles = 'velhrcurfile.nc'  ! Vel. cur. high freq. input file name 
    1261    velavadcpfiles = 'velavadcpfile.nc' ! Vel. ADCP daily av. input file name 
    1262    velhradcpfiles = 'velhradcpfile.nc' ! Vel. ADCP high freq. input file name 
    1263    velfbfiles = 'velfbfile.nc' ! Vel. feedback input observation file name 
    1264    dobsini = 20000101.000000  !  Initial date in window YYYYMMDD.HHMMSS 
    1265    dobsend = 20010101.000000  !  Final date in window YYYYMMDD.HHMMSS 
    1266    n1dint = 0  !               Type of vertical interpolation method 
    1267    n2dint = 0  !               Type of horizontal interpolation method 
    1268    ln_nea = .false.   !        Rejection of observations near land switch 
    1269    nmsshc     = 0     !        MSSH correction scheme 
    1270    mdtcorr = 1.61     !        MDT  correction 
    1271    mdtcutoff = 65.0   !        MDT cutoff for computed correction 
     1244                           !     ln_vel3d                Logical switch for velocity observations 
     1245                           !     ln_velavcur             Logical switch for velocity daily av. cur. 
     1246                           !     ln_velhrcur             Logical switch for velocity high freq. cur. 
     1247                           !     ln_velavadcp            Logical switch for velocity daily av. ADCP 
     1248                           !     ln_velhradcp            Logical switch for velocity high freq. ADCP 
     1249                           !     ln_velfb                Logical switch for feedback velocity data 
     1250                           !     ln_grid_global          Global distribtion of observations 
     1251                           !     ln_grid_search_lookup   Logical switch for obs grid search w/lookup table 
     1252                           !     grid_search_file        Grid search lookup file header 
     1253                           !     enactfiles              ENACT input observation file names 
     1254                           !     coriofiles              Coriolis input observation file name 
     1255   !                       ! profbfiles: Profile feedback input observation file name 
     1256   profbfiles = 'profiles_01.nc' 
     1257                           !     ln_profb_enatim         Enact feedback input time setting switch 
     1258                           !     slafilesact             Active SLA input observation file name 
     1259                           !     slafilespas             Passive SLA input observation file name 
     1260   !                       ! slafbfiles: Feedback SLA input observation file name 
     1261   slafbfiles = 'sla_01.nc' 
     1262                           !     sstfiles                GHRSST input observation file name 
     1263   !                       ! sstfbfiles: Feedback SST input observation file name 
     1264   sstfbfiles = 'sst_01.nc' 
     1265                           !     seaicefiles             Sea Ice input observation file names 
     1266   seaicefiles = 'seaice_01.nc' 
     1267                           !     velavcurfiles           Vel. cur. daily av. input file name 
     1268                           !     velhvcurfiles           Vel. cur. high freq. input file name 
     1269                           !     velavadcpfiles          Vel. ADCP daily av. input file name 
     1270                           !     velhvadcpfiles          Vel. ADCP high freq. input file name 
     1271                           !     velfbfiles              Vel. feedback input observation file name 
     1272                           !     dobsini                 Initial date in window YYYYMMDD.HHMMSS 
     1273                           !     dobsend                 Final date in window YYYYMMDD.HHMMSS 
     1274                           !     n1dint                  Type of vertical interpolation method 
     1275                           !     n2dint                  Type of horizontal interpolation method 
     1276                           !     ln_nea                  Rejection of observations near land switch 
     1277   nmsshc     = 0          ! MSSH correction scheme 
     1278                           !     mdtcorr                 MDT  correction 
     1279                           !     mdtcutoff               MDT cutoff for computed correction 
    12721280   ln_altbias = .false.    ! Logical switch for alt bias 
    12731281   ln_ignmis  = .true.     ! Logical switch for ignoring missing files 
    1274    endailyavtypes = 820    ! ENACT daily average types - array (use namelist_cfg to set more values) 
     1282                           !     endailyavtypes   ENACT daily average types 
    12751283   ln_grid_global = .true. 
    12761284   ln_grid_search_lookup = .false. 
    … …  
    12851293    ln_asmdin = .false.    !  Logical switch for Direct Initialization (DI) 
    12861294    ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU) 
     1295    ln_seaiceinc = .false. !  Logical switch for applying sea ice increments 
     1296    ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing 
    12871297    nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1] 
    12881298    nitdin    = 0          !  Timestep of background for DI in [0,nitend-nit000-1] 
    … …  
    13211331   rn_htrmax         =  200.0   ! max. depth of transition range 
    13221332/ 
     1333!----------------------------------------------------------------------- 
     1334&nambias   ! Bias pressure correctiom 
     1335!----------------------------------------------------------------------- 
     1336   ln_bias        = .false. 
     1337   ln_bias_asm    = .false. 
     1338   ln_bias_rlx    = .false. 
     1339   ln_bias_ofl    = .false. 
     1340   ln_bias_ts_app = .false. 
     1341   ln_bias_pc_app = .false.         
     1342   fb_t_asm       = 0.0 
     1343   fb_t_rlx       = 0.0 
     1344   fb_t_ofl       = 1.0 
     1345   fb_p_asm       = 1.0 
     1346   fb_p_rlx       = 1.0 
     1347   fb_p_ofl       = 0.0 
     1348   eft_rlx        = 365.0 
     1349   eft_asm        = 365.0 
     1350   t_rlx_upd      = 0.1 
     1351   t_asm_upd      = 0.1 
     1352   nn_lat_ramp    = 0           
     1353   bias_time_unit_asm = 86400.0 
     1354   bias_time_unit_rlx = 1.0 
     1355   bias_time_unit_ofl = 1.0  
     1356   cn_bias_tot    = "bias_tot.nc"  
     1357   cn_bias_asm    = "bias_asm.nc" 
     1358   cn_dir         = './'   
     1359   ln_bsyncro     = .FALSE.  
     1360   fctamp         = 1. 
     1361   rn_maxlat_bias = 23.0       
     1362   rn_minlat_bias = 10.0 
     1363   nn_bias_itwrt  = 15 
     1364   ln_itdecay     = .FALSE. 
     1365   ln_incpc       = .FALSE. 
     1366/ 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/CONFIG/SHARED/namelist_top_MEDUSA_et_al_ref

    r8280 r8657  
    6262   rn_ahtrc_0       =  2000.    !  horizontal eddy diffusivity for tracers [m2/s] 
    6363   rn_ahtrb_0       =     0.    !     background eddy diffusivity for ldf_iso [m2/s] 
     64   rn_fact_lap      =     1.    !     enhanced zonal eddy diffusivity 
    6465/ 
    6566!----------------------------------------------------------------------- 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r6486 r8657  
    119119            CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    120120#endif 
    121             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
     121!            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
     122            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
    122123            ! 
    123124            CALL iom_close( inum ) 
    … …  
    153154            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    154155            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     156            CALL iom_rstput( kt, nitdin_r, inum, 'avt'    , avt     ) 
    155157            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              ) 
    156158#if defined key_lim2 || defined key_lim3 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r7962 r8657  
    3939   USE ice_2            ! LIM2 
    4040#endif 
     41#if defined key_cice && defined key_asminc 
     42   USE sbc_ice, ONLY : & ! CICE Ice model variables 
     43   & ndaice_da, nfresh_da, nfsalt_da 
     44#endif 
    4145   USE sbc_oce          ! Surface boundary condition variables. 
    4246 
    … …  
    133137         &                 ln_asmdin, ln_asmiau,                           & 
    134138         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    135          &                 ln_salfix, salfixmin, nn_divdmp 
     139         &                 ln_salfix, salfixmin, nn_divdmp,                & 
     140         &                 ln_seaiceinc, ln_temnofreeze 
    136141      !!---------------------------------------------------------------------- 
    137142 
    … …  
    892897            ENDIF 
    893898 
     899         ELSE 
     900#if defined key_asminc 
     901            ssh_iau(:,:) = 0.0 
     902#endif 
    894903         ENDIF 
    895904 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r7747 r8657  
    8484      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace 
    8585      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace 
     86      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace 
    8687      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace 
    8788      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace 
    … …  
    9394      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse 
    9495      REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function 
    95       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace 
    9696 
    9797 
    … …  
    130130            zmask(:,:,:) = 0._wp 
    131131            zts(:,:,:,:) = 0._wp 
    132             zvn(:,:,:) = 0._wp 
    133132            DO jk = 1, jpkm1 
    134133               DO jj = 1, jpjm1 
    … …  
    138137                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid 
    139138                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc 
    140                      zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc 
    141139                  ENDDO 
    142140               ENDDO 
    … …  
    151149             tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1) 
    152150             sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1) 
    153              v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) ) 
     151             v_msf(:,:,1) = ptr_sjk( pvtr(:,:,:) ) 
    154152 
    155153             htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 ) 
    … …  
    177175                    tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    178176                    sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    179                     v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) )  
     177                    v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn) )  
    180178                    htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 ) 
    181179                    str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 ) 
    … …  
    202200             WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1) 
    203201             
    204             vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1)) 
     202            vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,1)) 
    205203            tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) ) 
    206204            tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) ) 
    … …  
    224222                    r1_sjk(:,1,jn) = 0._wp 
    225223                    WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn) 
    226                     vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn)) 
     224                    vsum = ptr_sj( pvtr(:,:,:), btmsk(:,:,jn)) 
    227225                    tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) 
    228226                    tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) 
    … …  
    408406            ENDIF 
    409407            IF( iom_use("zomsfeivglo") ) THEN 
    410                z3d(1,:,:) = ptr_sjk( v_eiv(:,:,:) )  ! zonal cumulative effective transport 
     408               DO jk=1,jpk 
     409                  DO jj=1,jpj 
     410                     DO ji=1,jpi 
     411                        zvn(ji,jj,jk) = v_eiv(ji,jj,jk) * fse3v(ji,jj,jk) * e1v(ji,jj) 
     412                     ENDDO 
     413                  ENDDO 
     414               ENDDO 
     415               z3d(1,:,:) = ptr_sjk( zvn(:,:,:) )  ! zonal cumulative effective transport 
    411416               DO jk = jpkm1,1,-1 
    412417                 z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)   ! effective j-Stream-Function (MSF) 
    … …  
    419424               IF( ln_subbas ) THEN 
    420425                  DO jn = 2, nptr                                    ! by sub-basins 
    421                      z3d(1,:,:) =  ptr_sjk( v_eiv(:,:,:), btmsk(:,:,jn) )  
     426                     z3d(1,:,:) =  ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) )  
    422427                     DO jk = jpkm1,1,-1 
    423428                        z3d(1,:,jk) = z3d(1,:,jk+1) - z3d(1,:,jk)    ! effective j-Stream-Function (MSF) 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8280 r8657  
    4747   USE iom 
    4848   USE ioipsl 
    49    USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities      
    50  
     49   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
     50   USE insitu_tem, ONLY: insitu_t, theta2t 
    5151#if defined key_lim2 
    5252   USE limwri_2  
    … …  
    164164       
    165165      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
     166      CALL theta2t ! in-situ temperature conversion 
     167      CALL iom_put( "tinsitu", insitu_t(:,:,:))    ! in-situ temperature 
    166168      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
    167169      IF ( iom_use("sbt") ) THEN 
    … …  
    202204         CALL iom_put( "taubot", z2d )            
    203205      ENDIF 
    204           
     206       
    205207      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
    206208      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r6486 r8657  
    355355         &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         & 
    356356         &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 
     357 
     358      ! Initilaise key variables at risk of being intercepted before properly set up.  
     359      e3t_0(:,:,:) = 0.0 
    357360         ! 
    358361#if defined key_vvl 
    … …  
    368371         &      ehu_b    (jpi,jpj)    , ehv_b  (jpi,jpj),                                                     & 
    369372         &      ehur_b   (jpi,jpj)    , ehvr_b (jpi,jpj),                                  STAT=ierr(5) )                           
     373 
     374      ! Initilaise key variables at risk of being intercepted before properly set up.  
     375      e3t_n(:,:,:) = 0.0 
    370376#endif 
    371377         ! 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r6486 r8657  
    4444   USE wrk_nemo        ! Memory Allocation 
    4545   USE timing          ! Timing 
     46   USE biaspar         ! bias correction variables 
    4647 
    4748   IMPLICIT NONE 
    … …  
    8485      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    8586      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     87      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st  ! tmp density storage for pressure corr 
     88      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_gru_st  ! tmp ua trends storage for pressure corr 
     89      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_grv_st  ! tmp va trends storage for pressure corr 
    8690      !!---------------------------------------------------------------------- 
    8791      ! 
    … …  
    9498      ENDIF 
    9599      ! 
     100      IF ( ln_bias .AND. ln_bias_pc_app ) THEN 
     101 
     102         !Allocate space for tempory variables 
     103         ALLOCATE( z_rhd_st(jpi,jpj,jpk), & 
     104            &      z_gru_st(jpi,jpj),     & 
     105            &      z_grv_st(jpi,jpj)      ) 
     106 
     107         z_rhd_st(:,:,:) = rhd(:,:,:)     ! store orig density  
     108         rhd(:,:,:)      = rhd_pc(:,:,:)  ! use pressure corrected density 
     109         z_gru_st(:,:)   = gru(:,:) 
     110         gru(:,:)        = gru_pc(:,:) 
     111         z_grv_st(:,:)   = grv(:,:) 
     112         grv(:,:)        = grv_pc(:,:) 
     113 
     114      ENDIF 
     115 
    96116      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    97117      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    … …  
    112132      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   & 
    113133         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     134      ! 
     135      IF ( ln_bias .AND. ln_bias_pc_app )  THEN 
     136         IF(lwp) THEN  
     137         WRITE(numout,*) " ! restore original density" 
     138         ENDIF 
     139         rhd(:,:,:) = z_rhd_st(:,:,:)     ! restore original density 
     140         gru(:,:)   = z_gru_st(:,:) 
     141         grv(:,:)   = z_grv_st(:,:) 
     142 
     143         !Deallocate tempory variables 
     144         DEALLOCATE( z_rhd_st,     & 
     145            &        z_gru_st,     & 
     146            &        z_grv_st      ) 
     147      ENDIF 
    114148      ! 
    115149      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg') 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r6487 r8657  
    7474      INTEGER, INTENT(in) ::   kt                      ! time step 
    7575      !  
    76       INTEGER             ::   jk                      ! dummy loop indice 
     76      INTEGER             ::   jk                      ! dummy loop indices 
    7777      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
    7878      !!---------------------------------------------------------------------- 
    … …  
    9494      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9595      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     96 
     97 
     98#if defined key_asminc 
     99      !                                                ! Include the IAU weighted SSH increment 
     100      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     101         CALL ssh_asm_inc( kt ) 
     102#if defined key_vvl 
     103! Don't directly adjust ssh but change hdivn at all levels instead 
     104! In trasbc also add in the heat and salt content associated with these changes at each level   
     105        DO jk = 1, jpkm1                                  
     106                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk)  
     107        END DO 
     108      ENDIF 
     109#endif 
     110#endif 
     111 
    96112 
    97113      !                                           !------------------------------! 
    … …  
    123139#endif 
    124140 
    125 #if defined key_asminc 
    126       !                                                ! Include the IAU weighted SSH increment 
    127       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    128          CALL ssh_asm_inc( kt ) 
    129          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 
    130       ENDIF 
    131 #endif 
    132141 
    133142      !                                           !------------------------------! 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r6498 r8657  
    110110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj 
    111111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_p, ht_p ! Meltpond fraction and depth 
     112    
     113   ! 
     114    
     115   ! 
     116#if defined key_asminc 
     117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ndaice_da          !: NEMO fresh water flux to ocean due to data assim 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfresh_da          !: NEMO salt flux to ocean due to data assim 
     119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfsalt_da          !: NEMO ice concentration change/second from data assim 
     120#endif 
     121       
    112122#endif 
    113123    
    … …  
    162172                ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , & 
    163173                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
     174#if defined key_asminc 
     175                ndaice_da(jpi,jpj)    , nfresh_da(jpi,jpj)    , nfsalt_da(jpi,jpj)    , & 
     176#endif 
    164177                sstfrz(jpi,jpj)       , STAT= ierr(1) ) 
    165178   ! Alex West: Allocating tn_ice with 5 categories.  When NEMO is used with CICE, this variable 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8280 r8657  
    21372137      REAL(wp) ::   zumax, zvmax 
    21382138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 
     2139      REAL(wp), POINTER, DIMENSION(:,:)   ::   zotx1_in, zoty1_in 
    21392140      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4    
    21402141      !!---------------------------------------------------------------------- 
    … …  
    21432144      ! 
    21442145      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 
     2146      CALL wrk_alloc( jpi,jpj, zotx1_in, zoty1_in) 
    21452147      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 ) 
    21462148 
    … …  
    24112413            zotx1(:,:) = un(:,:,1)   
    24122414            zoty1(:,:) = vn(:,:,1)   
    2413          ELSE         
     2415         ELSE 
    24142416            SELECT CASE( TRIM( sn_snd_crt%cldes ) ) 
    24152417            CASE( 'oce only'             )      ! C-grid ==> T 
    … …  
    25472549                  ENDDO 
    25482550               ENDDO 
    2549                 
     2551 
    25502552               ! Ensure any N fold and wrap columns are updated 
    25512553               CALL lbc_lnk(ztmp1, 'V', -1.0) 
    25522554               CALL lbc_lnk(ztmp2, 'U', -1.0) 
    2553                 
     2555                             
    25542556               ikchoix = -1 
    2555                CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix) 
     2557               ! We need copies of zotx1 and zoty2 in order to avoid problems  
     2558               ! caused by INTENTs used in the following subroutine.  
     2559               zotx1_in(:,:) = zotx1(:,:) 
     2560               zoty1_in(:,:) = zoty1(:,:) 
     2561               CALL repcmo (zotx1_in,ztmp2,ztmp1,zoty1_in,zotx1,zoty1,ikchoix) 
    25562562           ENDIF 
    25572563         ENDIF 
    … …  
    26222628      ! 
    26232629      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 
     2630      CALL wrk_dealloc( jpi,jpj, zotx1_in, zoty1_in ) 
    26242631      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 ) 
    26252632      ! 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r8280 r8657  
    5656                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             & 
    5757                flatn_f,fsurfn_f,fcondtopn_f,                    & 
     58#ifdef key_asminc 
     59                daice_da,fresh_da,fsalt_da,                    & 
     60#endif 
    5861                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
    5962                swvdr,swvdf,swidr,swidf,Tf,                      & 
    … …  
    301304  
    302305      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 ) 
     306 
     307#if defined key_asminc 
     308      ! Initialize fresh water and salt fluxes from data assim    
     309      !  and data assimilation index to cice  
     310      nfresh_da(:,:) = 0.0    
     311      nfsalt_da(:,:) = 0.0    
     312      ndaice_da(:,:) = 0.0          
     313#endif 
    303314      ! 
    304315      ! In coupled mode get extra fields from CICE for passing back to atmosphere 
    … …  
    454465      ENDIF 
    455466 
     467#if defined key_asminc 
     468!Ice concentration change (from assimilation) 
     469      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1) 
     470      Call nemo2cice(ztmp,daice_da,'T', 1. ) 
     471#endif  
     472 
    456473! Snowfall 
    457474! Ensure fsnow is positive (as in CICE routine prepare_forcing) 
    … …  
    716733         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt 
    717734      ENDIF 
     735 
     736#if defined key_asminc 
     737! Import fresh water and salt flux due to seaice da 
     738      CALL cice2nemo(fresh_da, nfresh_da,'T',1.0) 
     739      CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0) 
     740#endif 
    718741 
    719742! Release work space 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r7993 r8657  
    312312      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celcius] 
    313313      !                                                                ! 2 : salinity               [psu] 
    314       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
    315       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     314      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::   prd    ! in situ density            [-] 
     315      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::   prhop  ! potential density (surface referenced) 
    316316      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    317317      ! 
    … …  
    457457      END SELECT 
    458458      ! 
     459      CALL lbc_lnk( prd, 'T', 1.0_wp ) 
     460      ! 
    459461      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 
    460462      ! 
    … …  
    902904      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celcius,psu] 
    903905      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celcius-1,psu-1] 
    904       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
     906      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(inout) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
    905907      ! 
    906908      INTEGER  ::   ji, jj, jk      ! dummy loop indices 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r7771 r8657  
    549549      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    550550 
    551                                         !* sign of grad(H) at u- and v-points 
    552       mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
     551      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt 
     552      !!                                   !* sign of grad(H) at u- and v-points 
     553      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
     554      !! DO jj = 1, jpjm1 
     555      !!    DO ji = 1, jpim1 
     556      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     557      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     558      !!    END DO 
     559      !! END DO 
     560 
     561      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt 
     562                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0 
     563      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0 
    553564      DO jj = 1, jpjm1 
    554565         DO ji = 1, jpim1 
    555             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    556             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     566#if defined key_bbl_old_nonconserve 
     567             ! This key allows old (non conservative version) to be used for continuity of results 
     568             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     569             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     570#else 
     571            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
     572               mgrhu(ji,jj) = INT(  SIGN( 1.e0, & 
     573               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     574            ENDIF 
     575            !      
     576            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
     577               mgrhv(ji,jj) = INT(  SIGN( 1.e0, & 
     578               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     579            ENDIF 
     580#endif 
    557581         END DO 
    558582      END DO 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r7993 r8657  
    3333   USE timing          ! Timing 
    3434   USE eosbn2 
     35#if defined key_asminc    
     36   USE asminc          ! Assimilation increment 
     37#endif 
    3538 
    3639   IMPLICIT NONE 
    … …  
    120123      REAL(wp) ::   zfact, z1_e3t, zdep 
    121124      REAL(wp) ::   zalpha, zhk 
     125      REAL(wp) ::  zt_frz, zpress 
    122126      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    123127      !!---------------------------------------------------------------------- 
    … …  
    283287      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss 
    284288 
     289#if defined key_asminc 
     290! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM... 
     291! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0 
     292! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n) 
     293      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation 
     294         DO jj = 2, jpj  
     295            DO ji = fs_2, fs_jpim1 
     296               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) ) 
     297               DO jk = 1, jpkm1 
     298                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     299                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     300                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     301                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     302               END DO 
     303            END DO   
     304         END DO   
     305      ENDIF 
     306#endif 
     307  
    285308      IF( l_trdtra )   THEN                      ! send trends for further diagnostics 
    286309         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8280 r8657  
    227227#endif 
    228228      ! 
     229      ! Met Office addition: if failed, return non-zero exit code 
     230      IF( nstop /= 0 )  CALL exit( 9 )  
     231      ! 
    229232   END SUBROUTINE nemo_gcm 
    230233 
    … …  
    480483                            CALL dia_hsb_init   ! heat content, salt content and volume budgets 
    481484                            CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends 
     485                            CALL     bias_init  ! Pressure correction bias 
    482486      IF( lk_diaobs     ) THEN                  ! Observation & model comparison 
    483487                            CALL dia_obs_init            ! Initialize observational data 
    … …  
    646650      !!---------------------------------------------------------------------- 
    647651      USE diawri    , ONLY: dia_wri_alloc 
     652      USE insitu_tem, ONLY: insitu_tem_alloc 
    648653      USE dom_oce   , ONLY: dom_oce_alloc 
    649654      USE ldfdyn_oce, ONLY: ldfdyn_oce_alloc 
    … …  
    662667      ierr =        oce_alloc       ()          ! ocean 
    663668      ierr = ierr + dia_wri_alloc   () 
     669      ierr = ierr + insitu_tem_alloc() 
    664670      ierr = ierr + dom_oce_alloc   ()          ! ocean domain 
    665671      ierr = ierr + ldfdyn_oce_alloc()          ! ocean lateral  physics : dynamics 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8280 r8657  
    103103      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp 
    104104 
     105      IF( ln_bias )          CALL bias_opn( kstp ) 
     106 
    105107      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    106108      ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    … …  
    267269      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
    268270      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends 
     271      IF( ln_bias        )   CALL tra_bias   ( kstp ) 
    269272      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
    270273                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection 
    … …  
    290293               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    291294               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
     295            IF( ln_bias )    CALL dyn_bias( kstp ) 
    292296      ELSE                                                  ! centered hpg  (eos then time stepping) 
    293297         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    … …  
    303307         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    304308                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
     309         IF( ln_bias )       CALL dyn_bias( kstp ) 
    305310      ENDIF 
    306311 
    … …  
    377382      ENDIF 
    378383 
     384 
     385      IF( lrst_bias )          CALL bias_wrt     ( kstp ) 
     386 
    379387      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    380388      ! Coupled mode 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6491 r8657  
    100100 
    101101   USE crsfld           ! Standard output on coarse grid   (crs_fld routine) 
    102  
     102   USE biaspar          ! bias param 
     103   USE bias             ! bias routines                    (tra_bias routine 
     104                        !                                   dyn_bias routine) 
    103105   USE asminc           ! assimilation increments      (tra_asm_inc routine) 
    104106   !                                                   (dyn_asm_inc routine) 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/AGE/trcsms_age.F90

    r6715 r8657  
    5757      IF( nn_timing == 1 )  CALL timing_start('trc_sms_age') 
    5858      ! 
    59       IF(lwp) WRITE(numout,*) 
    60       IF(lwp) WRITE(numout,*) ' trc_sms_age:  AGE model' 
    61       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
     59      IF( kt == nittrc000 ) THEN 
     60         IF(lwp) WRITE(numout,*) 
     61         IF(lwp) WRITE(numout,*) ' trc_sms_age:  AGE model' 
     62         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
     63      ENDIF 
    6264 
    6365      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrage ) 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/C14b/trcnam_c14b.F90

    r6486 r8657  
    4949      ! definition of additional diagnostic as a structure 
    5050      INTEGER :: jl, jn 
    51       TYPE(DIAG), DIMENSION(jp_c14b_2d) :: c14dia2d 
    52       TYPE(DIAG), DIMENSION(jp_c14b_3d) :: c14dia3d 
    5351      !! 
    5452      NAMELIST/namc14date/ ndate_beg_b, nyear_res_b 
    55       NAMELIST/namc14dia/  c14dia2d, c14dia3d     ! additional diagnostics 
    5653      !!------------------------------------------------------------------- 
    5754      !                             ! Open namelist file 
    … …  
    7774      IF(lwp) WRITE(numout,*) '    initial year (aa)                  nyear_beg_b = ', nyear_beg_b 
    7875      ! 
    79       IF( .NOT.lk_iomput .AND. ln_diatrc ) THEN 
    80          ! 
    81          ! Namelist namc14dia 
    82          ! ------------------- 
    83          REWIND( numnatb_ref )              ! Namelist namc14dia in reference namelist : c14b diagnostics 
    84          READ  ( numnatb_ref, namc14dia, IOSTAT = ios, ERR = 903) 
    85 903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc14dia in reference namelist', lwp ) 
    86  
    87          REWIND( numnatb_cfg )              ! Namelist namc14dia in configuration namelist : c14b diagnostics 
    88          READ  ( numnatb_cfg, namc14dia, IOSTAT = ios, ERR = 904 ) 
    89 904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc14dia in configuration namelist', lwp ) 
    90          IF(lwm) WRITE ( numonb, namc14dia ) 
    91  
    92          DO jl = 1, jp_c14b_2d 
    93             jn = jp_c14b0_2d + jl - 1 
    94             ctrc2d(jn) = c14dia2d(jl)%sname 
    95             ctrc2l(jn) = c14dia2d(jl)%lname 
    96             ctrc2u(jn) = c14dia2d(jl)%units 
    97          END DO 
    98  
    99          DO jl = 1, jp_c14b_3d 
    100             jn = jp_c14b0_3d + jl - 1 
    101             ctrc3d(jn) = c14dia3d(jl)%sname 
    102             ctrc3l(jn) = c14dia3d(jl)%lname 
    103             ctrc3u(jn) = c14dia3d(jl)%units 
    104          END DO 
    105  
    106          IF(lwp) THEN                   ! control print 
    107             WRITE(numout,*) 
    108             WRITE(numout,*) ' Namelist : natadd' 
    109             DO jl = 1, jp_c14b_3d 
    110                jn = jp_c14b0_3d + jl - 1 
    111                WRITE(numout,*) '  3d diag nb : ', jn, '    short name : ', ctrc3d(jn), & 
    112                  &             '  long name  : ', ctrc3l(jn), '   unit : ', ctrc3u(jn) 
    113             END DO 
    114             WRITE(numout,*) ' ' 
    115  
    116             DO jl = 1, jp_c14b_2d 
    117                jn = jp_c14b0_2d + jl - 1 
    118                WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), & 
    119                  &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn) 
    120             END DO 
    121             WRITE(numout,*) ' ' 
    122          ENDIF 
    123          ! 
    124       ENDIF 
    12576 
    12677   IF(lwm) CALL FLUSH ( numonb )     ! flush output namelist C14b 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/CFC/trcnam_cfc.F90

    r8280 r8657  
    4747      INTEGER :: ios                 ! Local integer output status for namelist read 
    4848      INTEGER :: jl, jn 
    49       TYPE(DIAG), DIMENSION(jp_cfc_2d) :: cfcdia2d 
    5049      !! 
    5150      NAMELIST/namcfcdate/ ndate_beg, nyear_res, simu_type  
    52       NAMELIST/namcfcdia/  cfcdia2d     ! additional diagnostics 
    5351      !!---------------------------------------------------------------------- 
    5452      !                             ! Open namelist files 
    … …  
    8280      ! 
    8381 
    84       IF( .NOT.lk_iomput .AND. ln_diatrc ) THEN 
    85          ! 
    86          ! Namelist namcfcdia 
    87          ! ------------------- 
    88          REWIND( numnatc_ref )              ! Namelist namcfcdia in reference namelist : CFC diagnostics 
    89          READ  ( numnatc_ref, namcfcdia, IOSTAT = ios, ERR = 903) 
    90 903      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfcdia in reference namelist', lwp ) 
    91  
    92          REWIND( numnatc_cfg )              ! Namelist namcfcdia in configuration namelist : CFC diagnostics 
    93          READ  ( numnatc_cfg, namcfcdia, IOSTAT = ios, ERR = 904 ) 
    94 904      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfcdia in configuration namelist', lwp ) 
    95          IF(lwm) WRITE ( numonc, namcfcdia ) 
    96  
    97          DO jl = 1, jp_cfc_2d 
    98             jn = jp_cfc0_2d + jl - 1 
    99             ctrc2d(jn) = TRIM( cfcdia2d(jl)%sname ) 
    100             ctrc2l(jn) = TRIM( cfcdia2d(jl)%lname ) 
    101             ctrc2u(jn) = TRIM( cfcdia2d(jl)%units ) 
    102          END DO 
    103  
    104          IF(lwp) THEN                   ! control print 
    105             WRITE(numout,*) 
    106             WRITE(numout,*) ' Namelist : natadd' 
    107             DO jl = 1, jp_cfc_2d 
    108                jn = jp_cfc0_2d + jl - 1 
    109                WRITE(numout,*) '  2d diag nb : ', jn, '    short name : ', ctrc2d(jn), & 
    110                  &             '  long name  : ', ctrc2l(jn), '   unit : ', ctrc2u(jn) 
    111             END DO 
    112             WRITE(numout,*) ' ' 
    113          ENDIF 
    114          ! 
    115       ENDIF 
    116  
    11782   IF(lwm) CALL FLUSH ( numonc )     ! flush output namelist CFC 
    11883 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90

    r8280 r8657  
    257257      !ENDIF                                             
    258258      ! 
    259       IF( lk_iomput ) THEN 
    260          IF  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) ) 
    261          IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) ) 
    262          IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) ) 
    263          IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) ) 
    264          IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) ) 
    265          IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) ) 
    266       ELSE 
    267          IF( ln_diatrc ) THEN 
    268             trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1) 
    269             trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1) 
    270             trc2d(:,:,jp_cfc0_2d + 2) = qtr_cfc (:,:,2) 
    271             trc2d(:,:,jp_cfc0_2d + 3) = qint_cfc(:,:,2) 
    272             trc2d(:,:,jp_cfc0_2d + 4) = qtr_cfc (:,:,3) 
    273             trc2d(:,:,jp_cfc0_2d + 5) = qint_cfc(:,:,3) 
    274          END IF 
    275       END IF 
     259      IF  (iom_use("qtrCFC11"))  CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) ) 
     260      IF  (iom_use("qintCFC11")) CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) ) 
     261      IF  (iom_use("qtrCFC12"))  CALL iom_put( "qtrCFC12"  , qtr_cfc (:,:,2) ) 
     262      IF  (iom_use("qintCFC12")) CALL iom_put( "qintCFC12" , qint_cfc(:,:,2) ) 
     263      IF  (iom_use("qtrSF6"))    CALL iom_put( "qtrSF6"    , qtr_cfc (:,:,3) ) 
     264      IF  (iom_use("qintSF6"))   CALL iom_put( "qintSF6"   , qint_cfc(:,:,3) ) 
    276265      ! 
    277266      IF( l_trdtrc ) THEN 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/IDTRA/trcsms_idtra.F90

    r6829 r8657  
    165165      !ENDIF 
    166166      ! 
    167       IF( lk_iomput ) THEN 
    168167         CALL iom_put( "qtrIDTRA"  , qtr_idtra (:,:,1) ) 
    169168         CALL iom_put( "qintIDTRA" , qint_idtra(:,:,1) ) 
    170169         CALL iom_put( "invIDTRA" , inv_idtra(:,:,1) ) 
    171       ELSE 
    172          IF( ln_diatrc ) THEN 
    173             trc2d(:,:,jp_idtra0_2d    ) = qtr_idtra (:,:,1) 
    174             trc2d(:,:,jp_idtra0_2d + 1) = qint_idtra(:,:,1) 
    175             trc2d(:,:,jp_idtra0_2d + 2) = inv_idtra(:,:,1) 
    176          END IF 
    177       END IF 
    178170      ! 
    179171# if defined key_debug_medusa 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/MEDUSA/par_medusa.F90

    r6164 r8657  
    6363# endif 
    6464 
     65   ! assign an index in trc arrays for each PTS prognostic variables 
     66   INTEGER, PUBLIC, PARAMETER ::   jpchn_lc  =  1      !: non-diatom chlorophyll concentration 
     67   INTEGER, PUBLIC, PARAMETER ::   jpchd_lc  =  2      !: diatom     chlorophyll concentration 
     68   INTEGER, PUBLIC, PARAMETER ::   jpphn_lc  =  3      !: non-diatom concentration 
     69   INTEGER, PUBLIC, PARAMETER ::   jpphd_lc  =  4      !: diatom     concentration 
     70   INTEGER, PUBLIC, PARAMETER ::   jpzmi_lc  =  5      !: microzooplankton concentration 
     71   INTEGER, PUBLIC, PARAMETER ::   jpzme_lc  =  6      !: mesozooplankton  concentration 
     72   INTEGER, PUBLIC, PARAMETER ::   jpdin_lc  =  7      !: dissolved inorganic nitrogen concentration 
     73   INTEGER, PUBLIC, PARAMETER ::   jpsil_lc  =  8      !: silicic acid concentration 
     74   INTEGER, PUBLIC, PARAMETER ::   jpfer_lc  =  9      !: total iron concentration 
     75   INTEGER, PUBLIC, PARAMETER ::   jpdet_lc  =  10     !: slow-sinking detritus concentration 
     76   INTEGER, PUBLIC, PARAMETER ::   jppds_lc  =  11     !: diatom silicon concentration 
     77# if defined key_roam 
     78   INTEGER, PUBLIC, PARAMETER ::   jpdtc_lc  =  12     !: slow-sinking detritus carbon concentration 
     79   INTEGER, PUBLIC, PARAMETER ::   jpdic_lc  =  13     !: dissolved inorganic carbon concentration 
     80   INTEGER, PUBLIC, PARAMETER ::   jpalk_lc  =  14     !: alkalinity 
     81   INTEGER, PUBLIC, PARAMETER ::   jpoxy_lc  =  15     !: dissolved oxygen concentration 
     82# endif 
     83 
    6584#else 
    6685   !!--------------------------------------------------------------------- 

  • branches/NERC/dev_r5518_GO6_COAREbulk/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r8224 r8657  
    66   !! History : 
    77   !!  -   !  1999-07  (M. Levy)              original code 
    8    !!  -   !  2000-12  (E. Kestenare)         assign parameters to name individual tracers 
     8   !!  -   !  2000-12  (E. Kestenare)         assign parameters to name  
     9   !!                                         individual tracers 
    910   !!  -   !  2001-03  (M. Levy)              LNO3 + dia2d  
    1011   !! 2.0  !  2007-12  (C. Deltel, G. Madec)  F90 
    … …  
    1819   !!  -   !  2015-06  (A. Yool)              Update to include MOCSY 
    1920   !!  -   !  2015-07  (A. Yool)              Update for rolling averages 
    20    !!  -   !  2015-10  (J. Palm)              Update for diag outputs through iom_use   
     21   !!  -   !  2015-10  (J. Palm)              Update for diag outputs through  
     22   !!                                         iom_use   
    2123   !!  -   !  2016-11  (A. Yool)              Updated diags for CMIP6 
    2224   !!  -   !  2017-05  (A. Yool)              Added extra DMS calculation 
    … …  
    6062   !!   trc_bio_medusa        :   
    6163   !!---------------------------------------------------------------------- 
    62       USE oce_trc 
    63       USE trc 
    64       USE sms_medusa 
    65       USE lbclnk 
    66       USE prtctl_trc      ! Print control for debugging 
    67       USE trcsed_medusa 
    68       USE sbc_oce         ! surface forcing 
    69       USE sbcrnf          ! surface boundary condition: runoff variables 
    70       USE in_out_manager  ! I/O manager 
    71 # if defined key_iomput 
    72       USE iom 
    73       USE trcnam_medusa         ! JPALM 13-11-2015 -- if iom_use for diag 
    74       !!USE trc_nam_iom_medusa  ! JPALM 13-11-2015 -- if iom_use for diag 
    75 # endif 
    76 # if defined key_roam 
    77       USE gastransfer 
    78 #  if defined key_mocsy 
    79       USE mocsy_wrapper 
    80 #  else 
    81       USE trcco2_medusa 
    82 #  endif 
    83       USE trcoxy_medusa 
    84       !! Jpalm (08/08/2014) 
    85       USE trcdms_medusa 
    86 # endif 
    87       !! AXY (18/01/12): brought in for benthic timestepping 
    88       USE trcnam_trp      ! AXY (24/05/2013) 
    89       USE trdmxl_trc 
    90       USE trdtrc_oce  ! AXY (24/05/2013) 
    91  
    9264      !! AXY (30/01/14): necessary to find NaNs on HECTOR 
    9365      USE, INTRINSIC :: ieee_arithmetic  
    9466 
     67      USE bio_medusa_mod,             ONLY: b0, fdep1,                      &  
     68                                            ibenthic, idf, idfval,          & 
     69# if defined key_roam 
     70                                            f_xco2a,                        & 
     71                                            zalk, zdic, zoxy, zsal, ztmp,   & 
     72# endif 
     73# if defined key_mocsy 
     74                                            zpho,                           & 
     75# endif 
     76                                            zchd, zchn, zdet, zdin, zdtc,   & 
     77                                            zfer, zpds, zphd, zphn, zsil,   & 
     78                                            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 
     84      USE in_out_manager,             ONLY: lwp, numout, nn_date0 
     85# if defined key_iomput 
     86      USE iom,                        ONLY: lk_iomput 
     87# endif 
     88      USE lbclnk,                     ONLY: lbc_lnk 
     89      USE lib_mpp,                    ONLY: ctl_stop 
     90      USE oce,                        ONLY: tsb, tsn 
     91      USE par_kind,                   ONLY: wp 
     92      USE par_medusa,                 ONLY: jpalk, jpchd, jpchn, jpdet,     & 
     93                                            jpdic, jpdin, jpdtc, jpfer,     & 
     94                                            jpoxy, jppds, jpphd, jpphn,     & 
     95                                            jpsil, jpzme, jpzmi 
     96      USE par_oce,                    ONLY: jp_sal, jp_tem, jpi, jpim1,     & 
     97                                            jpj, jpjm1, jpk 
    9598      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    96       USE sbc_oce, ONLY: lk_oasis 
    97       USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, PCO2a_in_cpl, chloro_out_cpl 
     99      USE sbc_oce,                    ONLY: lk_oasis 
     100      USE sms_medusa,                 ONLY: hist_pco2 
     101      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
     102      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     103      USE carb_chem_mod,              ONLY: carb_chem 
     104      USE air_sea_mod,                ONLY: air_sea 
     105      USE plankton_mod,               ONLY: plankton 
     106      USE iron_chem_scav_mod,         ONLY: iron_chem_scav 
     107      USE detritus_mod,               ONLY: detritus 
     108      USE bio_medusa_update_mod,      ONLY: bio_medusa_update 
     109      USE bio_medusa_diag_mod,        ONLY: bio_medusa_diag 
     110      USE bio_medusa_diag_slice_mod,  ONLY: bio_medusa_diag_slice 
     111      USE bio_medusa_fin_mod,         ONLY: bio_medusa_fin 
    98112 
    99113      IMPLICIT NONE 
    100114      PRIVATE 
    101115       
    102       PUBLIC   trc_bio_medusa    ! called in ??? 
     116      PUBLIC   trc_bio_medusa    ! called in trcsms_medusa.F90 
    103117 
    104118   !!* Substitution 
    … …  
    113127 
    114128   SUBROUTINE trc_bio_medusa( kt ) 
    115       !!--------------------------------------------------------------------- 
     129      !!------------------------------------------------------------------ 
    116130      !!                     ***  ROUTINE trc_bio  *** 
    117131      !! 
    118       !! ** Purpose :   compute the now trend due to biogeochemical processes 
    119       !!              and add it to the general trend of passive tracers equations 
    120       !! 
    121       !! ** Method  :   each now biological flux is calculated in function of now 
    122       !!              concentrations of tracers. 
    123       !!              depending on the tracer, these fluxes are sources or sinks. 
    124       !!              the total of the sources and sinks for each tracer 
     132      !! ** Purpose : compute the now trend due to biogeochemical processes 
     133      !!              and add it to the general trend of passive tracers  
     134      !!              equations 
     135      !! 
     136      !! ** Method  : each now biological flux is calculated in function of 
     137      !!              now concentrations of tracers. 
     138      !!              depending on the tracer, these fluxes are sources or  
     139      !!              sinks. 
     140      !!              The total of the sources and sinks for each tracer 
    125141      !!              is added to the general trend. 
    126142      !!         
    … …  
    132148      !!              IF 'key_trc_diabio' defined , the biogeochemical trends 
    133149      !!              for passive tracers are saved for futher diagnostics. 
    134       !!--------------------------------------------------------------------- 
    135       !! 
    136       !! 
    137       !!----------------------------------------------------------------------             
     150      !!------------------------------------------------------------------ 
     151      !! 
     152      !! 
     153      !!------------------------------------------------------------------ 
    138154      !! Variable conventions 
    139       !!---------------------------------------------------------------------- 
     155      !!------------------------------------------------------------------ 
    140156      !! 
    141157      !! names: z*** - state variable  
    142       !!        f*** - function (or temporary variable used in part of a function) 
     158      !!        f*** - function (or temporary variable used in part of  
     159      !!               a function) 
    143160      !!        x*** - parameter 
    144161      !!        b*** - right-hand part (sources and sinks) 
    … …  
    151168      INTEGER  ::    ji,jj,jk,jn 
    152169      !! 
    153       !! AXY (27/07/10): add in indices for depth horizons (for sinking flux 
    154       !!                 and seafloor iron inputs) 
    155       !! INTEGER  ::    i0100, i0200, i0500, i1000, i1100 
    156       !! 
    157       !! model state variables 
    158       REAL(wp) ::    zchn,zchd,zphn,zphd,zpds,zzmi 
    159       REAL(wp) ::    zzme,zdet,zdtc,zdin,zsil,zfer 
    160       REAL(wp) ::    zage 
     170      INTEGER  ::    iball 
    161171# if defined key_roam 
    162       REAL(wp) ::    zdic, zalk, zoxy 
    163       REAL(wp) ::    ztmp, zsal 
    164 # endif 
    165 # if defined key_mocsy 
    166       REAL(wp) ::    zpho 
    167 # endif 
    168       !! 
    169       !! integrated source and sink terms 
    170       REAL(wp) ::    b0 
    171       !! AXY (23/08/13): changed from individual variables for each flux to 
    172       !!                 an array that holds all fluxes 
    173       REAL(wp), DIMENSION(jp_medusa) ::    btra 
    174       !! 
    175       !! primary production and chl related quantities       
    176       REAL(wp)                     ::    fthetan,faln,fchn1,fchn,fjln,fprn,frn 
    177       REAL(wp)                     ::    fthetad,fald,fchd1,fchd,fjld,fprd,frd 
    178       !! AXY (23/11/16): add in light-only limitation term (normalised 0-1 range) 
    179       REAL(wp)                     ::    fjlim_pn, fjlim_pd 
    180       !! AXY (03/02/11): add in Liebig terms 
    181       REAL(wp) ::    fpnlim, fpdlim 
    182       !! AXY (16/07/09): add in Eppley curve functionality 
    183       REAL(wp) ::    loc_T,fun_T,xvpnT,xvpdT 
    184       INTEGER  ::    ieppley 
    185       !! AXY (16/05/11): per Katya's prompting, add in new T-dependence 
    186       !!                 for phytoplankton growth only (i.e. no change 
    187       !!                 for remineralisation) 
    188       REAL(wp) ::    fun_Q10 
    189       !! AXY (01/03/10): add in mixed layer PP diagnostics 
    190       REAL(wp), DIMENSION(jpi,jpj) ::  fprn_ml,fprd_ml 
    191       !! 
    192       !! nutrient limiting factors 
    193       REAL(wp) ::    fnln,ffln            !! N and Fe 
    194       REAL(wp) ::    fnld,ffld,fsld,fsld2 !! N, Fe and Si 
    195       !! 
    196       !! silicon cycle 
    197       REAL(wp) ::    fsin,fnsi,fsin1,fnsi1,fnsi2,fprds,fsdiss 
    198       !! 
    199       !! iron cycle; includes parameters for Parekh et al. (2005) iron scheme 
    200       REAL(wp) ::    ffetop,ffebot,ffescav 
    201       REAL(wp) ::    xLgF, xFeT, xFeF, xFeL         !! state variables for iron-ligand system 
    202       REAL(wp), DIMENSION(jpi,jpj) ::  xFree        !! state variables for iron-ligand system 
    203       REAL(wp) ::    xb_coef_tmp, xb2M4ac           !! iron-ligand parameters 
    204       REAL(wp) ::    xmaxFeF,fdeltaFe               !! max Fe' parameters 
    205       !! 
    206       !! local parameters for Moore et al. (2004) alternative scavenging scheme 
    207       REAL(wp) ::    fbase_scav,fscal_sink,fscal_part,fscal_scav 
    208       !! 
    209       !! local parameters for Moore et al. (2008) alternative scavenging scheme 
    210       REAL(wp) ::    fscal_csink,fscal_sisink,fscal_casink 
    211       !! 
    212       !! local parameters for Galbraith et al. (2010) alternative scavenging scheme 
    213       REAL(wp) ::    xCscav1, xCscav2, xk_org, xORGscav  !! organic portion of scavenging 
    214       REAL(wp) ::    xk_inorg, xINORGscav                !! inorganic portion of scavenging 
    215       !! 
    216       !! microzooplankton grazing 
    217       REAL(wp) ::    fmi1,fmi,fgmipn,fgmid,fgmidc 
    218       REAL(wp) ::    finmi,ficmi,fstarmi,fmith,fmigrow,fmiexcr,fmiresp 
    219       !! 
    220       !! mesozooplankton grazing 
    221       REAL(wp) ::    fme1,fme,fgmepn,fgmepd,fgmepds,fgmezmi,fgmed,fgmedc 
    222       REAL(wp) ::    finme,ficme,fstarme,fmeth,fmegrow,fmeexcr,fmeresp 
    223       !! 
    224       !! mortality/Remineralisation (defunct parameter "fz" removed) 
    225       REAL(wp) ::    fdpn,fdpd,fdpds,fdzmi,fdzme,fdd 
    226 # if defined key_roam 
    227       REAL(wp) ::    fddc 
    228 # endif 
    229       REAL(wp) ::    fdpn2,fdpd2,fdpds2,fdzmi2,fdzme2 
    230       REAL(wp) ::    fslown, fslowc 
    231       REAL(wp), DIMENSION(jpi,jpj) ::    fslownflux, fslowcflux 
    232       REAL(wp) ::    fregen,fregensi 
    233       REAL(wp), DIMENSION(jpi,jpj) ::    fregenfast,fregenfastsi 
    234 # if defined key_roam 
    235       REAL(wp) ::    fregenc 
    236       REAL(wp), DIMENSION(jpi,jpj) ::    fregenfastc 
    237 # endif 
    238       !! 
    239       !! particle flux 
    240       REAL(WP) ::    fthk,fdep,fdep1,fdep2,flat,fcaco3 
    241       REAL(WP) ::    ftempn,ftempsi,ftempfe,ftempc,ftempca 
    242       REAL(wp) ::    freminn,freminsi,freminfe,freminc,freminca 
    243       REAL(wp), DIMENSION(jpi,jpj) ::    ffastn,ffastsi,ffastfe,ffastc,ffastca 
    244       REAL(wp) ::    fleftn,fleftsi,fleftfe,fleftc,fleftca 
    245       REAL(wp) ::    fheren,fheresi,fherefe,fherec,fhereca 
    246       REAL(wp) ::    fprotf 
    247       REAL(wp), DIMENSION(jpi,jpj) ::    fsedn,fsedsi,fsedfe,fsedc,fsedca 
    248       REAL(wp), DIMENSION(jpi,jpj) ::    fccd 
    249       REAL(wp) ::    fccd_dep 
    250       !! AXY (28/11/16): fix mbathy bug 
    251       INTEGER  ::    jmbathy 
    252       !! 
    253       !! AXY (06/07/11): alternative fast detritus schemes 
    254       REAL(wp) ::    fb_val, fl_sst 
    255       !! 
    256       !! AXY (08/07/11): fate of fast detritus reaching the seafloor 
    257       REAL(wp) ::    ffast2slown,ffast2slowfe,ffast2slowc 
    258       !! 
    259       !! conservation law 
    260       REAL(wp) ::    fnit0,fsil0,ffer0  
    261 # if defined key_roam 
    262       REAL(wp) ::    fcar0,falk0,foxy0  
    263 # endif       
     172      !! 
     173      INTEGER  ::    iyr1, iyr2 
     174      !! 
     175# endif 
    264176      !!  
    265177      !! temporary variables 
    266       REAL(wp) ::    fq0,fq1,fq2,fq3,fq4,fq5,fq6,fq7,fq8,fq9 
    267       !! 
    268       !! water column nutrient and flux integrals 
    269       REAL(wp), DIMENSION(jpi,jpj) ::    ftot_n,ftot_si,ftot_fe 
    270       REAL(wp), DIMENSION(jpi,jpj) ::    fflx_n,fflx_si,fflx_fe 
    271       REAL(wp), DIMENSION(jpi,jpj) ::    fifd_n,fifd_si,fifd_fe 
    272       REAL(wp), DIMENSION(jpi,jpj) ::    fofd_n,fofd_si,fofd_fe 
    273 # if defined key_roam 
    274       REAL(wp), DIMENSION(jpi,jpj) ::    ftot_c,ftot_a,ftot_o2 
    275       REAL(wp), DIMENSION(jpi,jpj) ::    fflx_c,fflx_a,fflx_o2 
    276       REAL(wp), DIMENSION(jpi,jpj) ::    fifd_c,fifd_a,fifd_o2 
    277       REAL(wp), DIMENSION(jpi,jpj) ::    fofd_c,fofd_a,fofd_o2 
    278 # endif 
    279       !! 
    280       !! zooplankton grazing integrals 
    281       REAL(wp), DIMENSION(jpi,jpj) ::    fzmi_i,fzmi_o,fzme_i,fzme_o 
    282       !! 
    283       !! limitation term temporary variables 
    284       REAL(wp), DIMENSION(jpi,jpj) ::    ftot_pn,ftot_pd 
    285       REAL(wp), DIMENSION(jpi,jpj) ::    ftot_zmi,ftot_zme,ftot_det,ftot_dtc 
    286       !! use ballast scheme (1) or simple exponential scheme (0; a conservation test) 
    287       INTEGER  ::    iball 
    288       !! use biological fluxes (1) or not (0) 
    289       INTEGER  ::    ibio_switch 
    290       !! 
    291       !! diagnose fluxes (should only be used in 1D runs) 
    292       INTEGER  ::    idf, idfval 
    293       !! 
    294       !! nitrogen and silicon production and consumption 
    295       REAL(wp) ::    fn_prod, fn_cons, fs_prod, fs_cons 
    296       REAL(wp), DIMENSION(jpi,jpj) ::    fnit_prod, fnit_cons, fsil_prod, fsil_cons 
    297 # if defined key_roam 
    298       !! 
    299       !! flags to help with calculating the position of the CCD 
    300       INTEGER, DIMENSION(jpi,jpj) ::     i2_omcal,i2_omarg 
    301       !! 
    302       !! ROAM air-sea flux and diagnostic parameters 
    303       REAL(wp) ::    f_wind 
    304       !! AXY (24/11/16): add xCO2 variable for atmosphere (what we actually have) 
    305       REAL(wp) ::    f_xco2a 
    306       REAL(wp) ::    f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_co2flux 
    307       REAL(wp) ::    f_TDIC, f_TALK, f_dcf, f_henry 
    308       REAL(wp) ::    f_uwind, f_vwind, f_pp0 
    309       REAL(wp) ::    f_kw660, f_o2flux, f_o2sat, f_o2sat3 
    310       REAL(wp), DIMENSION(jpi,jpj) ::    f_omcal, f_omarg 
    311       !! 
    312       !! AXY (23/06/15): additional diagnostics for MOCSY and oxygen 
    313       REAL(wp) ::    f_fco2w, f_BetaD, f_rhosw, f_opres, f_insitut, f_pco2atm, f_fco2atm 
    314       REAL(wp) ::    f_schmidtco2, f_kwco2, f_K0, f_co2starair, f_dpco2, f_kwo2 
    315       !! jpalm 14-07-2016: convert CO2flux diag from mmol/m2/d to kg/m2/s 
    316       REAL, PARAMETER :: weight_CO2_mol = 44.0095  !! g / mol 
    317       REAL, PARAMETER :: secs_in_day    = 86400.0  !! s / d 
    318       REAL, PARAMETER :: CO2flux_conv   = (1.e-6 * weight_CO2_mol) / secs_in_day 
    319  
    320       !! 
    321       INTEGER  ::    iters 
    322       REAL(wp) ::    f_year 
    323       INTEGER  ::    i_year 
    324       INTEGER  ::    iyr1, iyr2 
    325       !! 
    326       !! carbon, alkalinity production and consumption 
    327       REAL(wp) ::    fc_prod, fc_cons, fa_prod, fa_cons 
    328       REAL(wp), DIMENSION(jpi,jpj) ::    fcomm_resp 
    329       REAL(wp), DIMENSION(jpi,jpj) ::    fcar_prod, fcar_cons 
    330       !! 
    331       !! oxygen production and consumption (and non-consumption) 
    332       REAL(wp) ::    fo2_prod, fo2_cons, fo2_ncons, fo2_ccons 
    333       REAL(wp), DIMENSION(jpi,jpj) ::    foxy_prod, foxy_cons, foxy_anox 
    334       !! Jpalm (11-08-2014) 
    335       !! add DMS in MEDUSA for UKESM1 model 
    336       REAL(wp) ::    dms_surf 
    337       !! AXY (13/03/15): add in other DMS calculations 
    338       REAL(wp) ::    dms_andr, dms_simo, dms_aran, dms_hall, dms_andm, dms_nlim, dms_wtkn 
    339  
    340 # endif 
    341       !!  
    342       !! benthic fluxes 
    343       INTEGER  ::    ibenthic 
    344       REAL(wp), DIMENSION(jpi,jpj) :: f_sbenin_n, f_sbenin_fe,              f_sbenin_c 
    345       REAL(wp), DIMENSION(jpi,jpj) :: f_fbenin_n, f_fbenin_fe, f_fbenin_si, f_fbenin_c, f_fbenin_ca 
    346       REAL(wp), DIMENSION(jpi,jpj) :: f_benout_n, f_benout_fe, f_benout_si, f_benout_c, f_benout_ca 
    347       REAL(wp) ::    zfact 
    348       !!  
    349       !! benthic fluxes of CaCO3 that shouldn't happen because of lysocline 
    350       REAL(wp), DIMENSION(jpi,jpj) :: f_benout_lyso_ca 
    351       !! 
    352       !! riverine fluxes 
    353       REAL(wp), DIMENSION(jpi,jpj) :: f_runoff, f_riv_n, f_riv_si, f_riv_c, f_riv_alk 
    354       !! AXY (19/07/12): variables for local riverine fluxes to handle inputs below surface 
    355       REAL(wp) ::    f_riv_loc_n, f_riv_loc_si, f_riv_loc_c, f_riv_loc_alk 
    356       !! 
    357       !! Jpalm -- 11-10-2015 -- adapt diag to iom_use 
    358       !! 2D var for diagnostics. 
    359       REAL(wp), POINTER, DIMENSION(:,:  ) :: fprn2d, fdpn2d, fprd2d, fdpd2d, fprds2d, fsdiss2d, fgmipn2d 
    360       REAL(wp), POINTER, DIMENSION(:,:  ) :: fgmid2d, fdzmi2d, fgmepn2d, fgmepd2d, fgmezmi2d, fgmed2d 
    361       REAL(wp), POINTER, DIMENSION(:,:  ) :: fdzme2d, fslown2d, fdd2d, ffetop2d, ffebot2d, ffescav2d 
    362       REAL(wp), POINTER, DIMENSION(:,:  ) :: fjln2d, fnln2d, ffln2d, fjld2d, fnld2d, ffld2d, fsld2d2 
    363       REAL(wp), POINTER, DIMENSION(:,:  ) :: fsld2d, fregen2d, fregensi2d, ftempn2d, ftempsi2d, ftempfe2d 
    364       REAL(wp), POINTER, DIMENSION(:,:  ) :: ftempc2d, ftempca2d, freminn2d, freminsi2d, freminfe2d 
    365       REAL(wp), POINTER, DIMENSION(:,:  ) :: freminc2d, freminca2d 
    366       REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d 
    367 # if defined key_roam 
    368       REAL(wp), POINTER, DIMENSION(:,:  ) :: ffastca2d, rivn2d, rivsi2d, rivc2d, rivalk2d, fslowc2d 
    369       REAL(wp), POINTER, DIMENSION(:,:  ) :: fdpn22d, fdpd22d, fdzmi22d, fdzme22d, zimesn2d, zimesd2d 
    370       REAL(wp), POINTER, DIMENSION(:,:  ) :: zimesc2d, zimesdc2d, ziexcr2d, ziresp2d, zigrow2d, zemesn2d 
    371       REAL(wp), POINTER, DIMENSION(:,:  ) :: zemesd2d, zemesc2d, zemesdc2d, zeexcr2d, zeresp2d, zegrow2d 
    372       REAL(wp), POINTER, DIMENSION(:,:  ) :: mdetc2d, gmidc2d, gmedc2d, f_pco2a2d, f_pco2w2d, f_co2flux2d 
    373       REAL(wp), POINTER, DIMENSION(:,:  ) :: f_TDIC2d, f_TALK2d, f_kw6602d, f_pp02d, f_o2flux2d, f_o2sat2d 
    374       REAL(wp), POINTER, DIMENSION(:,:  ) :: dms_andr2d, dms_simo2d, dms_aran2d, dms_hall2d, dms_andm2d, dms_surf2d 
    375       REAL(wp), POINTER, DIMENSION(:,:  ) :: iben_n2d, iben_fe2d, iben_c2d, iben_si2d, iben_ca2d, oben_n2d 
    376       REAL(wp), POINTER, DIMENSION(:,:  ) :: oben_fe2d, oben_c2d, oben_si2d, oben_ca2d, sfr_ocal2d 
    377       REAL(wp), POINTER, DIMENSION(:,:  ) :: sfr_oarg2d, lyso_ca2d  
    378       !! AXY (23/11/16): extra MOCSY diagnostics 
    379       REAL(wp), POINTER, DIMENSION(:,:  ) :: f_xco2a_2d, f_fco2w_2d, f_fco2a_2d 
    380       REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnrhosw_2d, f_ocnschco2_2d, f_ocnkwco2_2d 
    381       REAL(wp), POINTER, DIMENSION(:,:  ) :: f_ocnk0_2d, f_co2starair_2d, f_ocndpco2_2d 
    382 # endif 
    383       !! 
    384       !! 3D var for diagnostics. 
    385       REAL(wp), POINTER, DIMENSION(:,:,:) :: tpp3d, detflux3d, remin3dn 
    386       !! 
    387 # if defined key_roam 
    388       !! AXY (04/11/16) 
    389       !! 2D var for new CMIP6 diagnostics (behind a key_roam ifdef for simplicity) 
    390       REAL(wp), POINTER, DIMENSION(:,:  ) :: fgco2, intdissic, intdissin, intdissisi, inttalk, o2min, zo2min 
    391       REAL(wp), POINTER, DIMENSION(:,:  ) :: fbddtalk, fbddtdic, fbddtdife, fbddtdin, fbddtdisi 
    392       !! 
    393       !! 3D var for new CMIP6 diagnostics 
    394       REAL(wp), POINTER, DIMENSION(:,:,:) :: tppd3 
    395       REAL(wp), POINTER, DIMENSION(:,:,:) :: bddtalk3, bddtdic3, bddtdife3, bddtdin3, bddtdisi3 
    396       REAL(wp), POINTER, DIMENSION(:,:,:) :: fd_nit3, fd_sil3, fd_car3, fd_cal3 
    397       REAL(wp), POINTER, DIMENSION(:,:,:) :: co33, co3satarag3, co3satcalc3, dcalc3 
    398       REAL(wp), POINTER, DIMENSION(:,:,:) :: expc3, expn3 
    399       REAL(wp), POINTER, DIMENSION(:,:,:) :: fediss3, fescav3 
    400       REAL(wp), POINTER, DIMENSION(:,:,:) :: migrazp3, migrazd3, megrazp3, megrazd3, megrazz3 
    401       REAL(wp), POINTER, DIMENSION(:,:,:) :: o2sat3, pbsi3, pcal3, remoc3 
    402       REAL(wp), POINTER, DIMENSION(:,:,:) :: pnlimj3, pnlimn3, pnlimfe3, pdlimj3, pdlimn3, pdlimfe3, pdlimsi3 
    403 # endif 
    404       !!--------------------------------------------------------------------- 
     178      REAL(wp) ::    fq0,fq1,fq2,fq3,fq4 
     179      !! 
     180      !!------------------------------------------------------------------ 
    405181 
    406182# if defined key_debug_medusa 
    … …  
    421197      ibenthic = 1 
    422198 
    423       !! not sure what this is for; it's not used anywhere; commenting out 
    424       !! fbodn(:,:) = 0.e0    
    425  
    426       !! 
    427       IF( ln_diatrc ) THEN 
    428          !! blank 2D diagnostic array 
    429          trc2d(:,:,:) = 0.e0 
    430          !! 
    431          !! blank 3D diagnostic array 
    432          trc3d(:,:,:,:) = 0.e0 
    433       ENDIF 
    434  
    435       !!---------------------------------------------------------------------- 
     199      !!------------------------------------------------------------------ 
    436200      !! b0 is present for debugging purposes; using b0 = 0 sets the tendency 
    437201      !! terms of all biological equations to 0. 
    438       !!---------------------------------------------------------------------- 
     202      !!------------------------------------------------------------------ 
    439203      !! 
    440204      !! AXY (03/09/14): probably not the smartest move ever, but it'll fit 
    … …  
    446210      b0 = 1. 
    447211# endif 
    448       !!---------------------------------------------------------------------- 
     212      !!------------------------------------------------------------------ 
    449213      !! fast detritus ballast scheme (0 = no; 1 = yes) 
    450214      !! alternative to ballast scheme is same scheme but with no ballast 
    451215      !! protection (not dissimilar to Martin et al., 1987) 
    452       !!---------------------------------------------------------------------- 
     216      !!------------------------------------------------------------------ 
    453217      !! 
    454218      iball = 1 
    455219 
    456       !!---------------------------------------------------------------------- 
     220      !!------------------------------------------------------------------ 
    457221      !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output 
    458222      !! these should *only* be used in 1D since they give comprehensive 
    459223      !! output for ecological functions in the model; primarily used in 
    460224      !! debugging 
    461       !!---------------------------------------------------------------------- 
     225      !!------------------------------------------------------------------ 
    462226      !! 
    463227      idf    = 0 
    … …  
    470234      endif 
    471235 
    472       !!---------------------------------------------------------------------- 
    473       !! blank fast-sinking detritus 2D fields 
    474       !!---------------------------------------------------------------------- 
    475       !! 
    476       ffastn(:,:)  = 0.0        !! organic nitrogen 
    477       ffastsi(:,:) = 0.0        !! biogenic silicon 
    478       ffastfe(:,:) = 0.0        !! organic iron 
    479       ffastc(:,:)  = 0.0        !! organic carbon 
    480       ffastca(:,:) = 0.0        !! biogenic calcium carbonate 
    481       !! 
    482       fsedn(:,:)   = 0.0        !! Seafloor flux of N  
    483       fsedsi(:,:)  = 0.0        !! Seafloor flux of Si 
    484       fsedfe(:,:)  = 0.0        !! Seafloor flux of Fe 
    485       fsedc(:,:)   = 0.0        !! Seafloor flux of C 
    486       fsedca(:,:)  = 0.0        !! Seafloor flux of CaCO3 
    487       !! 
    488       fregenfast(:,:)   = 0.0   !! integrated  N regeneration (fast detritus) 
    489       fregenfastsi(:,:) = 0.0   !! integrated Si regeneration (fast detritus) 
    490 # if defined key_roam 
    491       fregenfastc(:,:)  = 0.0   !! integrated  C regeneration (fast detritus) 
    492 # endif 
    493       !! 
    494       fccd(:,:)    = 0.0        !! last depth level before CCD 
    495  
    496       !!---------------------------------------------------------------------- 
    497       !! blank nutrient/flux inventories 
    498       !!---------------------------------------------------------------------- 
    499       !! 
    500       fflx_n(:,:)  = 0.0        !! nitrogen flux total 
    501       fflx_si(:,:) = 0.0        !! silicon  flux total 
    502       fflx_fe(:,:) = 0.0        !! iron     flux total 
    503       fifd_n(:,:)  = 0.0        !! nitrogen fast detritus production 
    504       fifd_si(:,:) = 0.0        !! silicon  fast detritus production 
    505       fifd_fe(:,:) = 0.0        !! iron     fast detritus production 
    506       fofd_n(:,:)  = 0.0        !! nitrogen fast detritus remineralisation 
    507       fofd_si(:,:) = 0.0        !! silicon  fast detritus remineralisation 
    508       fofd_fe(:,:) = 0.0        !! iron     fast detritus remineralisation 
    509 # if defined key_roam 
    510       fflx_c(:,:)  = 0.0        !! carbon     flux total 
    511       fflx_a(:,:)  = 0.0        !! alkalinity flux total 
    512       fflx_o2(:,:) = 0.0        !! oxygen     flux total 
    513       ftot_c(:,:)  = 0.0        !! carbon     inventory 
    514       ftot_a(:,:)  = 0.0        !! alkalinity inventory 
    515       ftot_o2(:,:) = 0.0        !! oxygen     inventory 
    516       fifd_c(:,:)  = 0.0        !! carbon     fast detritus production 
    517       fifd_a(:,:)  = 0.0        !! alkalinity fast detritus production 
    518       fifd_o2(:,:) = 0.0        !! oxygen     fast detritus production 
    519       fofd_c(:,:)  = 0.0        !! carbon     fast detritus remineralisation 
    520       fofd_a(:,:)  = 0.0        !! alkalinity fast detritus remineralisation 
    521       fofd_o2(:,:) = 0.0        !! oxygen     fast detritus remineralisation 
    522       !! 
    523       fnit_prod(:,:) = 0.0      !! (organic)   nitrogen production 
    524       fnit_cons(:,:) = 0.0      !! (organic)   nitrogen consumption 
    525       fsil_prod(:,:) = 0.0      !! (inorganic) silicon production 
    526       fsil_cons(:,:) = 0.0      !! (inorganic) silicon consumption 
    527       fcar_prod(:,:) = 0.0      !! (organic)   carbon production 
    528       fcar_cons(:,:) = 0.0      !! (organic)   carbon consumption 
    529       !! 
    530       foxy_prod(:,:) = 0.0      !! oxygen production 
    531       foxy_cons(:,:) = 0.0      !! oxygen consumption 
    532       foxy_anox(:,:) = 0.0      !! unrealised oxygen consumption 
    533       !! 
    534 # endif 
    535       ftot_n(:,:)   = 0.0       !! N inventory  
    536       ftot_si(:,:)  = 0.0       !! Si inventory 
    537       ftot_fe(:,:)  = 0.0       !! Fe inventory 
    538       ftot_pn(:,:)  = 0.0       !! integrated non-diatom phytoplankton 
    539       ftot_pd(:,:)  = 0.0       !! integrated diatom     phytoplankton 
    540       ftot_zmi(:,:) = 0.0       !! integrated microzooplankton 
    541       ftot_zme(:,:) = 0.0       !! integrated mesozooplankton 
    542       ftot_det(:,:) = 0.0       !! integrated slow detritus, nitrogen 
    543       ftot_dtc(:,:) = 0.0       !! integrated slow detritus, carbon 
    544       !! 
    545       fzmi_i(:,:)  = 0.0        !! material grazed by microzooplankton 
    546       fzmi_o(:,:)  = 0.0        !! ... sum of fate of this material 
    547       fzme_i(:,:)  = 0.0        !! material grazed by  mesozooplankton 
    548       fzme_o(:,:)  = 0.0        !! ... sum of fate of this material 
    549       !! 
    550       f_sbenin_n(:,:)  = 0.0    !! slow detritus N  -> benthic pool 
    551       f_sbenin_fe(:,:) = 0.0    !! slow detritus Fe -> benthic pool 
    552       f_sbenin_c(:,:)  = 0.0    !! slow detritus C  -> benthic pool 
    553       f_fbenin_n(:,:)  = 0.0    !! fast detritus N  -> benthic pool 
    554       f_fbenin_fe(:,:) = 0.0    !! fast detritus Fe -> benthic pool 
    555       f_fbenin_si(:,:) = 0.0    !! fast detritus Si -> benthic pool 
    556       f_fbenin_c(:,:)  = 0.0    !! fast detritus C  -> benthic pool 
    557       f_fbenin_ca(:,:) = 0.0    !! fast detritus Ca -> benthic pool 
    558       !! 
    559       f_benout_n(:,:)  = 0.0    !! benthic N  pool  -> dissolved 
    560       f_benout_fe(:,:) = 0.0    !! benthic Fe pool  -> dissolved 
    561       f_benout_si(:,:) = 0.0    !! benthic Si pool  -> dissolved 
    562       f_benout_c(:,:)  = 0.0    !! benthic C  pool  -> dissolved 
    563       f_benout_ca(:,:) = 0.0    !! benthic Ca pool  -> dissolved 
    564       !! 
    565       f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool  -> dissolved (when it shouldn't!) 
    566       !! 
    567       f_runoff(:,:)  = 0.0      !! riverine runoff 
    568       f_riv_n(:,:)   = 0.0      !! riverine N   input  
    569       f_riv_si(:,:)  = 0.0      !! riverine Si  input  
    570       f_riv_c(:,:)   = 0.0      !! riverine C   input  
    571       f_riv_alk(:,:) = 0.0      !! riverine alk input  
    572       !!  
    573       !! Jpalm -- 06-03-2017 -- Forgotten var to init 
    574       f_omarg(:,:) = 0.0        !! 
    575       f_omcal(:,:) = 0.0  
    576       xFree(:,:) = 0.0          !! state variables for iron-ligand system 
    577       fcomm_resp(:,:) = 0.0  
    578       fprn_ml(:,:) = 0.0        !! mixed layer PP diagnostics 
    579       fprd_ml(:,:) = 0.0        !! mixed layer PP diagnostics 
    580       !! 
    581       fslownflux(:,:) = 0.0 
    582       fslowcflux(:,:) = 0.0 
    583  
    584       !! 
    585       !! allocate and initiate 2D diag 
    586       !! ----------------------------- 
    587       !! Juju :: add kt condition !! 
    588       IF ( lk_iomput .AND. .NOT.  ln_diatrc ) THEN  
    589          !! 
    590          if ( kt == nittrc000 )   CALL trc_nam_iom_medusa !! initialise iom_use test 
    591          !! 
    592          CALL wrk_alloc( jpi, jpj,      zw2d ) 
    593          zw2d(:,:)      = 0.0   !! 
    594          IF ( med_diag%PRN%dgsave ) THEN 
    595             CALL wrk_alloc( jpi, jpj,   fprn2d    ) 
    596             fprn2d(:,:)      = 0.0 !! 
    597          ENDIF 
    598          IF ( med_diag%MPN%dgsave ) THEN 
    599             CALL wrk_alloc( jpi, jpj,   fdpn2d    ) 
    600             fdpn2d(:,:)      = 0.0 !! 
    601          ENDIF 
    602          IF ( med_diag%PRD%dgsave ) THEN 
    603             CALL wrk_alloc( jpi, jpj,   fprd2d    ) 
    604             fprd2d(:,:)      = 0.0 !! 
    605          ENDIF 
    606          IF( med_diag%MPD%dgsave ) THEN 
    607             CALL wrk_alloc( jpi, jpj,   fdpd2d    ) 
    608             fdpd2d(:,:)      = 0.0 !! 
    609          ENDIF 
    610          IF( med_diag%OPAL%dgsave ) THEN 
    611             CALL wrk_alloc( jpi, jpj,   fprds2d    ) 
    612             fprds2d(:,:)      = 0.0 !! 
    613          ENDIF 
    614          IF( med_diag%OPALDISS%dgsave ) THEN 
    615             CALL wrk_alloc( jpi, jpj,   fsdiss2d    ) 
    616             fsdiss2d(:,:)      = 0.0 !! 
    617          ENDIF 
    618          IF( med_diag%GMIPn%dgsave ) THEN 
    619             CALL wrk_alloc( jpi, jpj,   fgmipn2d    ) 
    620             fgmipn2d(:,:)      = 0.0 !! 
    621          ENDIF 
    622          IF( med_diag%GMID%dgsave ) THEN 
    623             CALL wrk_alloc( jpi, jpj,   fgmid2d    ) 
    624             fgmid2d(:,:)      = 0.0 !! 
    625          ENDIF 
    626          IF( med_diag%MZMI%dgsave ) THEN 
    627             CALL wrk_alloc( jpi, jpj,   fdzmi2d    ) 
    628             fdzmi2d(:,:)      = 0.0 !! 
    629          ENDIF 
    630          IF( med_diag%GMEPN%dgsave ) THEN 
    631             CALL wrk_alloc( jpi, jpj,   fgmepn2d    ) 
    632             fgmepn2d(:,:)      = 0.0 !! 
    633          ENDIF 
    634          IF( med_diag%GMEPD%dgsave ) THEN 
    635             CALL wrk_alloc( jpi, jpj,   fgmepd2d    ) 
    636             fgmepd2d(:,:)      = 0.0 !! 
    637          ENDIF 
    638          IF( med_diag%GMEZMI%dgsave ) THEN 
    639             CALL wrk_alloc( jpi, jpj,   fgmezmi2d    ) 
    640             fgmezmi2d(:,:)      = 0.0 !! 
    641          ENDIF 
    642          IF( med_diag%GMED%dgsave ) THEN 
    643             CALL wrk_alloc( jpi, jpj,   fgmed2d    ) 
    644             fgmed2d(:,:)      = 0.0 !! 
    645          ENDIF 
    646          IF( med_diag%MZME%dgsave ) THEN 
    647             CALL wrk_alloc( jpi, jpj,   fdzme2d    ) 
    648             fdzme2d(:,:)      = 0.0 !! 
    649          ENDIF 
    650          IF( med_diag%DETN%dgsave ) THEN 
    651             CALL wrk_alloc( jpi, jpj,   fslown2d    ) 
    652             fslown2d(:,:)      = 0.0 !! 
    653          ENDIF 
    654          IF( med_diag%MDET%dgsave ) THEN 
    655             CALL wrk_alloc( jpi, jpj,   fdd2d    ) 
    656             fdd2d(:,:)      = 0.0 !! 
    657          ENDIF       
    658          IF( med_diag%AEOLIAN%dgsave ) THEN 
    659             CALL wrk_alloc( jpi, jpj,   ffetop2d    ) 
    660             ffetop2d(:,:)      = 0.0 !! 
    661          ENDIF 
    662          IF( med_diag%BENTHIC%dgsave ) THEN 
    663             CALL wrk_alloc( jpi, jpj,    ffebot2d   ) 
    664             ffebot2d(:,:)      = 0.0 !! 
    665          ENDIF 
    666          IF( med_diag%SCAVENGE%dgsave ) THEN 
    667             CALL wrk_alloc( jpi, jpj,   ffescav2d    ) 
    668             ffescav2d(:,:)      = 0.0 !! 
    669          ENDIF 
    670          IF( med_diag%PN_JLIM%dgsave ) THEN 
    671             CALL wrk_alloc( jpi, jpj,   fjln2d    ) 
    672             fjln2d(:,:)      = 0.0 !! 
    673          ENDIF 
    674          IF( med_diag%PN_NLIM%dgsave ) THEN 
    675             CALL wrk_alloc( jpi, jpj,   fnln2d    ) 
    676             fnln2d(:,:)      = 0.0 !! 
    677          ENDIF 
    678          IF( med_diag%PN_FELIM%dgsave ) THEN 
    679             CALL wrk_alloc( jpi, jpj,   ffln2d    ) 
    680             ffln2d(:,:)      = 0.0 !! 
    681          ENDIF 
    682          IF( med_diag%PD_JLIM%dgsave ) THEN 
    683             CALL wrk_alloc( jpi, jpj,   fjld2d    ) 
    684             fjld2d(:,:)      = 0.0 !! 
    685          ENDIF 
    686          IF( med_diag%PD_NLIM%dgsave ) THEN 
    687             CALL wrk_alloc( jpi, jpj,   fnld2d    ) 
    688             fnld2d(:,:)      = 0.0 !! 
    689          ENDIF 
    690          IF( med_diag%PD_FELIM%dgsave ) THEN 
    691             CALL wrk_alloc( jpi, jpj,   ffld2d    ) 
    692             ffld2d(:,:)      = 0.0 !! 
    693          ENDIF 
    694          IF( med_diag%PD_SILIM%dgsave ) THEN 
    695             CALL wrk_alloc( jpi, jpj,   fsld2d2    ) 
    696             fsld2d2(:,:)      = 0.0 !! 
    697          ENDIF 
    698          IF( med_diag%PDSILIM2%dgsave ) THEN 
    699             CALL wrk_alloc( jpi, jpj,   fsld2d    ) 
    700             fsld2d(:,:)      = 0.0 !! 
    701          ENDIF 
    702 !! 
    703 !! skip SDT_XXXX diagnostics here 
    704 !! 
    705          IF( med_diag%TOTREG_N%dgsave ) THEN 
    706             CALL wrk_alloc( jpi, jpj,   fregen2d    ) 
    707             fregen2d(:,:)      = 0.0 !! 
    708          ENDIF 
    709          IF( med_diag%TOTRG_SI%dgsave ) THEN 
    710             CALL wrk_alloc( jpi, jpj,   fregensi2d    ) 
    711             fregensi2d(:,:)      = 0.0 !! 
    712          ENDIF 
    713 !! 
    714 !! skip REG_XXXX diagnostics here 
    715 !! 
    716          IF( med_diag%FASTN%dgsave ) THEN 
    717             CALL wrk_alloc( jpi, jpj,   ftempn2d    ) 
    718             ftempn2d(:,:)      = 0.0 !! 
    719          ENDIF 
    720          IF( med_diag%FASTSI%dgsave ) THEN 
    721             CALL wrk_alloc( jpi, jpj,   ftempsi2d    ) 
    722             ftempsi2d(:,:)      = 0.0 !! 
    723          ENDIF 
    724          IF( med_diag%FASTFE%dgsave ) THEN 
    725             CALL wrk_alloc( jpi, jpj,  ftempfe2d     ) 
    726             ftempfe2d(:,:)      = 0.0 !! 
    727          ENDIF 
    728          IF( med_diag%FASTC%dgsave ) THEN 
    729             CALL wrk_alloc( jpi, jpj,  ftempc2d     ) 
    730             ftempc2d(:,:)      = 0.0 !! 
    731          ENDIF 
    732          IF( med_diag%FASTCA%dgsave ) THEN 
    733             CALL wrk_alloc( jpi, jpj,   ftempca2d    ) 
    734             ftempca2d(:,:)      = 0.0 !! 
    735          ENDIF      
    736 !! 
    737 !! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here 
    738 !! 
    739          IF( med_diag%REMINN%dgsave ) THEN 
    740             CALL wrk_alloc( jpi, jpj,    freminn2d   ) 
    741             freminn2d(:,:)      = 0.0 !! 
    742          ENDIF 
    743          IF( med_diag%REMINSI%dgsave ) THEN 
    744             CALL wrk_alloc( jpi, jpj,    freminsi2d   ) 
    745             freminsi2d(:,:)      = 0.0 !! 
    746          ENDIF 
    747          IF( med_diag%REMINFE%dgsave ) THEN 
    748             CALL wrk_alloc( jpi, jpj,    freminfe2d   ) 
    749             freminfe2d(:,:)      = 0.0 !! 
    750          ENDIF 
    751          IF( med_diag%REMINC%dgsave ) THEN 
    752             CALL wrk_alloc( jpi, jpj,   freminc2d    ) 
    753             freminc2d(:,:)      = 0.0 !!  
    754          ENDIF 
    755          IF( med_diag%REMINCA%dgsave ) THEN 
    756             CALL wrk_alloc( jpi, jpj,   freminca2d    ) 
    757             freminca2d(:,:)      = 0.0 !! 
    758          ENDIF 
    759 # if defined key_roam 
    760 !! 
    761 !! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here 
    762 !! 
    763          IF( med_diag%RR_0100%dgsave ) THEN 
    764             CALL wrk_alloc( jpi, jpj,    ffastca2d   ) 
    765             ffastca2d(:,:)      = 0.0 !! 
    766          ENDIF 
    767  
    768          IF( med_diag%ATM_PCO2%dgsave ) THEN 
    769             CALL wrk_alloc( jpi, jpj,    f_pco2a2d   ) 
    770             f_pco2a2d(:,:)      = 0.0 !! 
    771          ENDIF 
    772 !! 
    773 !! skip OCN_PH diagnostic here 
    774 !! 
    775          IF( med_diag%OCN_PCO2%dgsave ) THEN 
    776             CALL wrk_alloc( jpi, jpj,    f_pco2w2d   ) 
    777             f_pco2w2d(:,:)      = 0.0 !! 
    778          ENDIF 
    779 !! 
    780 !! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here 
    781 !! 
    782          IF( med_diag%CO2FLUX%dgsave ) THEN 
    783             CALL wrk_alloc( jpi, jpj,   f_co2flux2d    ) 
    784             f_co2flux2d(:,:)      = 0.0 !! 
    785          ENDIF 
    786 !! 
    787 !! skip OM_XXX diagnostics here 
    788 !! 
    789          IF( med_diag%TCO2%dgsave ) THEN 
    790             CALL wrk_alloc( jpi, jpj,   f_TDIC2d    ) 
    791             f_TDIC2d(:,:)      = 0.0 !! 
    792          ENDIF 
    793          IF( med_diag%TALK%dgsave ) THEN 
    794             CALL wrk_alloc( jpi, jpj,    f_TALK2d   ) 
    795             f_TALK2d(:,:)      = 0.0 !! 
    796          ENDIF 
    797          IF( med_diag%KW660%dgsave ) THEN 
    798             CALL wrk_alloc( jpi, jpj,    f_kw6602d   ) 
    799             f_kw6602d(:,:)      = 0.0 !! 
    800          ENDIF 
    801          IF( med_diag%ATM_PP0%dgsave ) THEN 
    802             CALL wrk_alloc( jpi, jpj,    f_pp02d   ) 
    803             f_pp02d(:,:)      = 0.0 !! 
    804          ENDIF 
    805          IF( med_diag%O2FLUX%dgsave ) THEN 
    806             CALL wrk_alloc( jpi, jpj,   f_o2flux2d    ) 
    807             f_o2flux2d(:,:)      = 0.0 !! 
    808          ENDIF 
    809          IF( med_diag%O2SAT%dgsave ) THEN 
    810             CALL wrk_alloc( jpi, jpj,    f_o2sat2d   ) 
    811             f_o2sat2d(:,:)      = 0.0 !! 
    812          ENDIF  
    813 !! 
    814 !! skip XXX_CCD diagnostics here 
    815 !!  
    816          IF( med_diag%SFR_OCAL%dgsave ) THEN 
    817             CALL wrk_alloc( jpi, jpj,    sfr_ocal2d  ) 
    818             sfr_ocal2d(:,:)      = 0.0 !! 
    819          ENDIF 
    820          IF( med_diag%SFR_OARG%dgsave ) THEN 
    821             CALL wrk_alloc( jpi, jpj,    sfr_oarg2d  ) 
    822             sfr_oarg2d(:,:)      = 0.0 !! 
    823          ENDIF 
    824 !! 
    825 !! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here 
    826 !!  
    827          IF( med_diag%IBEN_N%dgsave ) THEN 
    828             CALL wrk_alloc( jpi, jpj,    iben_n2d  ) 
    829             iben_n2d(:,:)      = 0.0 !! 
    830          ENDIF 
    831          IF( med_diag%IBEN_FE%dgsave ) THEN 
    832             CALL wrk_alloc( jpi, jpj,   iben_fe2d   ) 
    833             iben_fe2d(:,:)      = 0.0 !! 
    834          ENDIF 
    835          IF( med_diag%IBEN_C%dgsave ) THEN 
    836             CALL wrk_alloc( jpi, jpj,   iben_c2d   ) 
    837             iben_c2d(:,:)      = 0.0 !! 
    838          ENDIF 
    839          IF( med_diag%IBEN_SI%dgsave ) THEN 
    840             CALL wrk_alloc( jpi, jpj,   iben_si2d   ) 
    841             iben_si2d(:,:)      = 0.0 !! 
    842          ENDIF 
    843          IF( med_diag%IBEN_CA%dgsave ) THEN 
    844             CALL wrk_alloc( jpi, jpj,   iben_ca2d   ) 
    845             iben_ca2d(:,:)      = 0.0 !! 
    846          ENDIF 
    847          IF( med_diag%OBEN_N%dgsave ) THEN 
    848             CALL wrk_alloc( jpi, jpj,    oben_n2d  ) 
    849             oben_n2d(:,:)      = 0.0 !! 
    850          ENDIF 
    851          IF( med_diag%OBEN_FE%dgsave ) THEN 
    852             CALL wrk_alloc( jpi, jpj,    oben_fe2d  ) 
    853             oben_fe2d(:,:)      = 0.0 !! 
    854          ENDIF 
    855          IF( med_diag%OBEN_C%dgsave ) THEN 
    856             CALL wrk_alloc( jpi, jpj,    oben_c2d  ) 
    857             oben_c2d(:,:)      = 0.0 !! 
    858          ENDIF 
    859          IF( med_diag%OBEN_SI%dgsave ) THEN 
    860             CALL wrk_alloc( jpi, jpj,    oben_si2d  ) 
    861             oben_si2d(:,:)      = 0.0 !! 
    862          ENDIF 
    863          IF( med_diag%OBEN_CA%dgsave ) THEN 
    864             CALL wrk_alloc( jpi, jpj,    oben_ca2d  ) 
    865             oben_ca2d(:,:)      = 0.0 !! 
    866          ENDIF 
    867 !! 
    868 !! skip BEN_XX diagnostics here 
    869 !! 
    870          IF( med_diag%RIV_N%dgsave ) THEN 
    871             CALL wrk_alloc( jpi, jpj,    rivn2d   ) 
    872             rivn2d(:,:)      = 0.0 !! 
    873          ENDIF 
    874          IF( med_diag%RIV_SI%dgsave ) THEN 
    875             CALL wrk_alloc( jpi, jpj,    rivsi2d   ) 
    876             rivsi2d(:,:)      = 0.0 !! 
    877          ENDIF 
    878          IF( med_diag%RIV_C%dgsave ) THEN 
    879             CALL wrk_alloc( jpi, jpj,   rivc2d    ) 
    880             rivc2d(:,:)      = 0.0 !! 
    881          ENDIF 
    882          IF( med_diag%RIV_ALK%dgsave ) THEN 
    883             CALL wrk_alloc( jpi, jpj,    rivalk2d   ) 
    884             rivalk2d(:,:)      = 0.0 !! 
    885          ENDIF 
    886          IF( med_diag%DETC%dgsave ) THEN 
    887             CALL wrk_alloc( jpi, jpj,    fslowc2d   ) 
    888             fslowc2d(:,:)      = 0.0 !! 
    889          ENDIF  
    890 !! 
    891 !! skip SDC_XXXX, INVTXXX diagnostics here 
    892 !! 
    893          IF( med_diag%LYSO_CA%dgsave ) THEN 
    894             CALL wrk_alloc( jpi, jpj,    lyso_ca2d  ) 
    895             lyso_ca2d(:,:)      = 0.0 !! 
    896          ENDIF 
    897 !! 
    898 !! skip COM_RESP diagnostic here 
    899 !! 
    900          IF( med_diag%PN_LLOSS%dgsave ) THEN 
    901             CALL wrk_alloc( jpi, jpj,    fdpn22d   ) 
    902             fdpn22d(:,:)      = 0.0 !! 
    903          ENDIF 
    904          IF( med_diag%PD_LLOSS%dgsave ) THEN 
    905             CALL wrk_alloc( jpi, jpj,    fdpd22d   ) 
    906             fdpd22d(:,:)      = 0.0 !! 
    907          ENDIF 
    908          IF( med_diag%ZI_LLOSS%dgsave ) THEN 
    909             CALL wrk_alloc( jpi, jpj,    fdzmi22d   ) 
    910             fdzmi22d(:,:)      = 0.0 !! 
    911          ENDIF 
    912          IF( med_diag%ZE_LLOSS%dgsave ) THEN 
    913             CALL wrk_alloc( jpi, jpj,   fdzme22d    ) 
    914             fdzme22d(:,:)      = 0.0 !! 
    915          ENDIF 
    916          IF( med_diag%ZI_MES_N%dgsave ) THEN    
    917             CALL wrk_alloc( jpi, jpj,   zimesn2d    ) 
    918             zimesn2d(:,:)      = 0.0 !! 
    919          ENDIF 
    920          IF( med_diag%ZI_MES_D%dgsave ) THEN 
    921             CALL wrk_alloc( jpi, jpj,    zimesd2d   ) 
    922             zimesd2d(:,:)      = 0.0 !! 
    923          ENDIF 
    924          IF( med_diag%ZI_MES_C%dgsave ) THEN 
    925             CALL wrk_alloc( jpi, jpj,    zimesc2d   ) 
    926             zimesc2d(:,:)      = 0.0 !! 
    927          ENDIF 
    928          IF( med_diag%ZI_MESDC%dgsave ) THEN 
    929             CALL wrk_alloc( jpi, jpj,    zimesdc2d   ) 
    930             zimesdc2d(:,:)      = 0.0 !! 
    931          ENDIF 
    932          IF( med_diag%ZI_EXCR%dgsave ) THEN 
    933             CALL wrk_alloc( jpi, jpj,     ziexcr2d  ) 
    934             ziexcr2d(:,:)      = 0.0 !! 
    935          ENDIF 
    936          IF( med_diag%ZI_RESP%dgsave ) THEN 
    937             CALL wrk_alloc( jpi, jpj,    ziresp2d   ) 
    938             ziresp2d(:,:)      = 0.0 !! 
    939          ENDIF 
    940          IF( med_diag%ZI_GROW%dgsave ) THEN 
    941             CALL wrk_alloc( jpi, jpj,    zigrow2d   ) 
    942             zigrow2d(:,:)      = 0.0 !! 
    943          ENDIF 
    944          IF( med_diag%ZE_MES_N%dgsave ) THEN 
    945             CALL wrk_alloc( jpi, jpj,   zemesn2d    ) 
    946             zemesn2d(:,:)      = 0.0 !! 
    947          ENDIF 
    948          IF( med_diag%ZE_MES_D%dgsave ) THEN 
    949             CALL wrk_alloc( jpi, jpj,    zemesd2d   ) 
    950             zemesd2d(:,:)      = 0.0 !! 
    951          ENDIF 
    952          IF( med_diag%ZE_MES_C%dgsave ) THEN 
    953             CALL wrk_alloc( jpi, jpj,    zemesc2d   ) 
    954             zemesc2d(:,:)      = 0.0 !! 
    955          ENDIF 
    956          IF( med_diag%ZE_MESDC%dgsave ) THEN 
    957             CALL wrk_alloc( jpi, jpj,    zemesdc2d   ) 
    958             zemesdc2d(:,:)      = 0.0 !! 
    959          ENDIF 
    960          IF( med_diag%ZE_EXCR%dgsave ) THEN 
    961             CALL wrk_alloc( jpi, jpj,    zeexcr2d   ) 
    962             zeexcr2d(:,:)      = 0.0 !! 
    963          ENDIF                   
    964          IF( med_diag%ZE_RESP%dgsave ) THEN 
    965             CALL wrk_alloc( jpi, jpj,    zeresp2d   ) 
    966             zeresp2d(:,:)      = 0.0 !! 
    967          ENDIF 
    968          IF( med_diag%ZE_GROW%dgsave ) THEN 
    969             CALL wrk_alloc( jpi, jpj,    zegrow2d   ) 
    970             zegrow2d(:,:)      = 0.0 !! 
    971          ENDIF 
    972          IF( med_diag%MDETC%dgsave ) THEN 
    973             CALL wrk_alloc( jpi, jpj,   mdetc2d    ) 
    974             mdetc2d(:,:)      = 0.0 !! 
    975          ENDIF 
    976          IF( med_diag%GMIDC%dgsave ) THEN 
    977             CALL wrk_alloc( jpi, jpj,    gmidc2d   ) 
    978             gmidc2d(:,:)      = 0.0 !! 
    979          ENDIF 
    980          IF( med_diag%GMEDC%dgsave ) THEN 
    981             CALL wrk_alloc( jpi, jpj,    gmedc2d   ) 
    982             gmedc2d(:,:)      = 0.0 !! 
    983          ENDIF 
    984 !! 
    985 !! skip INT_XXX diagnostics here 
    986 !! 
    987          IF (jdms .eq. 1) THEN 
    988             IF( med_diag%DMS_SURF%dgsave ) THEN 
    989                CALL wrk_alloc( jpi, jpj,   dms_surf2d    ) 
    990                dms_surf2d(:,:)      = 0.0 !! 
    991             ENDIF 
    992             IF( med_diag%DMS_ANDR%dgsave ) THEN 
    993                CALL wrk_alloc( jpi, jpj,   dms_andr2d    ) 
    994                dms_andr2d(:,:)      = 0.0 !! 
    995             ENDIF 
    996             IF( med_diag%DMS_SIMO%dgsave ) THEN 
    997                CALL wrk_alloc( jpi, jpj,  dms_simo2d     ) 
    998                dms_simo2d(:,:)      = 0.0 !! 
    999             ENDIF 
    1000             IF( med_diag%DMS_ARAN%dgsave ) THEN 
    1001                CALL wrk_alloc( jpi, jpj,   dms_aran2d    ) 
    1002                dms_aran2d(:,:)      = 0.0 !! 
    1003             ENDIF 
    1004             IF( med_diag%DMS_HALL%dgsave ) THEN 
    1005                CALL wrk_alloc( jpi, jpj,   dms_hall2d    ) 
    1006                dms_hall2d(:,:)      = 0.0 !! 
    1007             ENDIF 
    1008             IF( med_diag%DMS_ANDM%dgsave ) THEN 
    1009                CALL wrk_alloc( jpi, jpj,   dms_andm2d    ) 
    1010                dms_andm2d(:,:)      = 0.0 !! 
    1011             ENDIF 
    1012          ENDIF    
    1013          !! 
    1014          !! AXY (24/11/16): extra MOCSY diagnostics, 2D 
    1015          IF( med_diag%ATM_XCO2%dgsave ) THEN 
    1016             CALL wrk_alloc( jpi, jpj, f_xco2a_2d      ) 
    1017             f_xco2a_2d(:,:)      = 0.0 !! 
    1018          ENDIF 
    1019          IF( med_diag%OCN_FCO2%dgsave ) THEN 
    1020             CALL wrk_alloc( jpi, jpj, f_fco2w_2d      ) 
    1021             f_fco2w_2d(:,:)      = 0.0 !! 
    1022          ENDIF 
    1023          IF( med_diag%ATM_FCO2%dgsave ) THEN 
    1024             CALL wrk_alloc( jpi, jpj, f_fco2a_2d      ) 
    1025             f_fco2a_2d(:,:)      = 0.0 !! 
    1026          ENDIF 
    1027          IF( med_diag%OCN_RHOSW%dgsave ) THEN 
    1028             CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d   ) 
    1029             f_ocnrhosw_2d(:,:)      = 0.0 !! 
    1030          ENDIF 
    1031          IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
    1032             CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d  ) 
    1033             f_ocnschco2_2d(:,:)      = 0.0 !! 
    1034          ENDIF 
    1035          IF( med_diag%OCN_KWCO2%dgsave ) THEN 
    1036             CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d   ) 
    1037             f_ocnkwco2_2d(:,:)      = 0.0 !! 
    1038          ENDIF 
    1039          IF( med_diag%OCN_K0%dgsave ) THEN 
    1040             CALL wrk_alloc( jpi, jpj, f_ocnk0_2d      ) 
    1041             f_ocnk0_2d(:,:)      = 0.0 !! 
    1042          ENDIF 
    1043          IF( med_diag%CO2STARAIR%dgsave ) THEN 
    1044             CALL wrk_alloc( jpi, jpj, f_co2starair_2d ) 
    1045             f_co2starair_2d(:,:)      = 0.0 !! 
    1046          ENDIF 
    1047          IF( med_diag%OCN_DPCO2%dgsave ) THEN 
    1048             CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d   ) 
    1049             f_ocndpco2_2d(:,:)      = 0.0 !! 
    1050          ENDIF 
    1051 # endif   
    1052          IF( med_diag%TPP3%dgsave ) THEN 
    1053             CALL wrk_alloc( jpi, jpj, jpk,       tpp3d ) 
    1054             tpp3d(:,:,:)      = 0.0 !!  
    1055          ENDIF 
    1056          IF( med_diag%DETFLUX3%dgsave ) THEN 
    1057             CALL wrk_alloc( jpi, jpj, jpk,        detflux3d ) 
    1058             detflux3d(:,:,:)      = 0.0 !!  
    1059          ENDIF 
    1060          IF( med_diag%REMIN3N%dgsave ) THEN 
    1061              CALL wrk_alloc( jpi, jpj, jpk,        remin3dn ) 
    1062              remin3dn(:,:,:)      = 0.0 !!  
    1063           ENDIF 
    1064           !!  
    1065           !! AXY (10/11/16): CMIP6 diagnostics, 2D 
    1066           !! JPALM -- 17-11-16 -- put fgco2 alloc out of diag request 
    1067           !!                   needed for coupling/passed through restart 
    1068           !! IF( med_diag%FGCO2%dgsave ) THEN 
    1069              CALL wrk_alloc( jpi, jpj,   fgco2    ) 
    1070              fgco2(:,:)      = 0.0 !! 
    1071           !! ENDIF 
    1072           IF( med_diag%INTDISSIC%dgsave ) THEN 
    1073              CALL wrk_alloc( jpi, jpj,   intdissic    ) 
    1074              intdissic(:,:)  = 0.0 !! 
    1075           ENDIF           
    1076           IF( med_diag%INTDISSIN%dgsave ) THEN 
    1077              CALL wrk_alloc( jpi, jpj,   intdissin    ) 
    1078              intdissin(:,:)  = 0.0 !! 
    1079           ENDIF           
    1080           IF( med_diag%INTDISSISI%dgsave ) THEN 
    1081              CALL wrk_alloc( jpi, jpj,   intdissisi    ) 
    1082              intdissisi(:,:)  = 0.0 !! 
    1083           ENDIF           
    1084           IF( med_diag%INTTALK%dgsave ) THEN 
    1085              CALL wrk_alloc( jpi, jpj,   inttalk    ) 
    1086              inttalk(:,:)  = 0.0 !! 
    1087           ENDIF           
    1088           IF( med_diag%O2min%dgsave ) THEN 
    1089              CALL wrk_alloc( jpi, jpj,   o2min    ) 
    1090              o2min(:,:)  = 1.e3 !! set to high value as we're looking for min(o2) 
    1091           ENDIF           
    1092           IF( med_diag%ZO2min%dgsave ) THEN 
    1093              CALL wrk_alloc( jpi, jpj,   zo2min    ) 
    1094              zo2min(:,:)  = 0.0 !! 
    1095           ENDIF           
    1096           IF( med_diag%FBDDTALK%dgsave  ) THEN 
    1097              CALL wrk_alloc( jpi, jpj, fbddtalk  ) 
    1098              fbddtalk(:,:)  = 0.0 !!  
    1099           ENDIF 
    1100           IF( med_diag%FBDDTDIC%dgsave  ) THEN 
    1101              CALL wrk_alloc( jpi, jpj, fbddtdic  ) 
    1102              fbddtdic(:,:)  = 0.0 !!  
    1103           ENDIF 
    1104           IF( med_diag%FBDDTDIFE%dgsave ) THEN 
    1105              CALL wrk_alloc( jpi, jpj, fbddtdife ) 
    1106              fbddtdife(:,:) = 0.0 !!  
    1107           ENDIF 
    1108           IF( med_diag%FBDDTDIN%dgsave  ) THEN 
    1109              CALL wrk_alloc( jpi, jpj, fbddtdin  ) 
    1110              fbddtdin(:,:)  = 0.0 !!  
    1111           ENDIF 
    1112           IF( med_diag%FBDDTDISI%dgsave ) THEN 
    1113              CALL wrk_alloc( jpi, jpj, fbddtdisi ) 
    1114              fbddtdisi(:,:) = 0.0 !!  
    1115           ENDIF 
    1116           !!  
    1117           !! AXY (10/11/16): CMIP6 diagnostics, 3D 
    1118           IF( med_diag%TPPD3%dgsave     ) THEN 
    1119              CALL wrk_alloc( jpi, jpj, jpk, tppd3     ) 
    1120              tppd3(:,:,:)     = 0.0 !!  
    1121           ENDIF 
    1122           IF( med_diag%BDDTALK3%dgsave  ) THEN 
    1123              CALL wrk_alloc( jpi, jpj, jpk, bddtalk3  ) 
    1124              bddtalk3(:,:,:)  = 0.0 !!  
    1125           ENDIF 
    1126           IF( med_diag%BDDTDIC3%dgsave  ) THEN 
    1127              CALL wrk_alloc( jpi, jpj, jpk, bddtdic3  ) 
    1128              bddtdic3(:,:,:)  = 0.0 !!  
    1129           ENDIF 
    1130           IF( med_diag%BDDTDIFE3%dgsave ) THEN 
    1131              CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 ) 
    1132              bddtdife3(:,:,:) = 0.0 !!  
    1133           ENDIF 
    1134           IF( med_diag%BDDTDIN3%dgsave  ) THEN 
    1135              CALL wrk_alloc( jpi, jpj, jpk, bddtdin3  ) 
    1136              bddtdin3(:,:,:)  = 0.0 !!  
    1137           ENDIF 
    1138           IF( med_diag%BDDTDISI3%dgsave ) THEN 
    1139              CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 ) 
    1140              bddtdisi3(:,:,:) = 0.0 !!  
    1141           ENDIF 
    1142           IF( med_diag%FD_NIT3%dgsave   ) THEN 
    1143              CALL wrk_alloc( jpi, jpj, jpk, fd_nit3   ) 
    1144              fd_nit3(:,:,:)   = 0.0 !!  
    1145           ENDIF 
    1146           IF( med_diag%FD_SIL3%dgsave   ) THEN 
    1147              CALL wrk_alloc( jpi, jpj, jpk, fd_sil3   ) 
    1148              fd_sil3(:,:,:)   = 0.0 !!  
    1149           ENDIF 
    1150           IF( med_diag%FD_CAR3%dgsave   ) THEN 
    1151              CALL wrk_alloc( jpi, jpj, jpk, fd_car3   ) 
    1152              fd_car3(:,:,:)   = 0.0 !!  
    1153           ENDIF 
    1154           IF( med_diag%FD_CAL3%dgsave   ) THEN 
    1155              CALL wrk_alloc( jpi, jpj, jpk, fd_cal3   ) 
    1156              fd_cal3(:,:,:)   = 0.0 !!  
    1157           ENDIF 
    1158           IF( med_diag%DCALC3%dgsave    ) THEN 
    1159              CALL wrk_alloc( jpi, jpj, jpk, dcalc3    ) 
    1160              dcalc3(:,:,: )   = 0.0 !!  
    1161           ENDIF 
    1162           IF( med_diag%EXPC3%dgsave     ) THEN 
    1163              CALL wrk_alloc( jpi, jpj, jpk, expc3   ) 
    1164              expc3(:,:,: )    = 0.0 !!  
    1165           ENDIF 
    1166           IF( med_diag%EXPN3%dgsave     ) THEN 
    1167              CALL wrk_alloc( jpi, jpj, jpk, expn3   ) 
    1168              expn3(:,:,: )    = 0.0 !!  
    1169           ENDIF 
    1170           IF( med_diag%FEDISS3%dgsave   ) THEN 
    1171              CALL wrk_alloc( jpi, jpj, jpk, fediss3   ) 
    1172              fediss3(:,:,: )  = 0.0 !!  
    1173           ENDIF 
    1174           IF( med_diag%FESCAV3%dgsave   ) THEN 
    1175              CALL wrk_alloc( jpi, jpj, jpk, fescav3   ) 
    1176              fescav3(:,:,: )  = 0.0 !!  
    1177           ENDIF 
    1178           IF( med_diag%MIGRAZP3%dgsave   ) THEN 
    1179              CALL wrk_alloc( jpi, jpj, jpk, migrazp3  ) 
    1180              migrazp3(:,:,: )  = 0.0 !!  
    1181           ENDIF 
    1182           IF( med_diag%MIGRAZD3%dgsave   ) THEN 
    1183              CALL wrk_alloc( jpi, jpj, jpk, migrazd3  ) 
    1184              migrazd3(:,:,: )  = 0.0 !!  
    1185           ENDIF 
    1186           IF( med_diag%MEGRAZP3%dgsave   ) THEN 
    1187              CALL wrk_alloc( jpi, jpj, jpk, megrazp3  ) 
    1188              megrazp3(:,:,: )  = 0.0 !!  
    1189           ENDIF 
    1190           IF( med_diag%MEGRAZD3%dgsave   ) THEN 
    1191              CALL wrk_alloc( jpi, jpj, jpk, megrazd3  ) 
    1192              megrazd3(:,:,: )  = 0.0 !!  
    1193           ENDIF 
    1194           IF( med_diag%MEGRAZZ3%dgsave   ) THEN 
    1195              CALL wrk_alloc( jpi, jpj, jpk, megrazz3  ) 
    1196              megrazz3(:,:,: )  = 0.0 !!  
    1197           ENDIF 
    1198           IF( med_diag%O2SAT3%dgsave     ) THEN 
    1199              CALL wrk_alloc( jpi, jpj, jpk, o2sat3    ) 
    1200              o2sat3(:,:,: )    = 0.0 !!  
    1201           ENDIF 
    1202           IF( med_diag%PBSI3%dgsave      ) THEN 
    1203              CALL wrk_alloc( jpi, jpj, jpk, pbsi3     ) 
    1204              pbsi3(:,:,: )     = 0.0 !!  
    1205           ENDIF 
    1206           IF( med_diag%PCAL3%dgsave      ) THEN 
    1207              CALL wrk_alloc( jpi, jpj, jpk, pcal3     ) 
    1208              pcal3(:,:,: )     = 0.0 !!  
    1209           ENDIF 
    1210           IF( med_diag%REMOC3%dgsave     ) THEN 
    1211              CALL wrk_alloc( jpi, jpj, jpk, remoc3    ) 
    1212              remoc3(:,:,: )    = 0.0 !!  
    1213           ENDIF 
    1214           IF( med_diag%PNLIMJ3%dgsave    ) THEN 
    1215              CALL wrk_alloc( jpi, jpj, jpk, pnlimj3   ) 
    1216              pnlimj3(:,:,: )   = 0.0 !!  
    1217           ENDIF 
    1218           IF( med_diag%PNLIMN3%dgsave    ) THEN 
    1219              CALL wrk_alloc( jpi, jpj, jpk, pnlimn3   ) 
    1220              pnlimn3(:,:,: )   = 0.0 !!  
    1221           ENDIF 
    1222           IF( med_diag%PNLIMFE3%dgsave   ) THEN 
    1223              CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3  ) 
    1224              pnlimfe3(:,:,: )  = 0.0 !!  
    1225           ENDIF 
    1226           IF( med_diag%PDLIMJ3%dgsave    ) THEN 
    1227              CALL wrk_alloc( jpi, jpj, jpk, pdlimj3   ) 
    1228              pdlimj3(:,:,: )   = 0.0 !!  
    1229           ENDIF 
    1230           IF( med_diag%PDLIMN3%dgsave    ) THEN 
    1231              CALL wrk_alloc( jpi, jpj, jpk, pdlimn3   ) 
    1232              pdlimn3(:,:,: )   = 0.0 !!  
    1233           ENDIF 
    1234           IF( med_diag%PDLIMFE3%dgsave   ) THEN 
    1235              CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3  ) 
    1236              pdlimfe3(:,:,: )  = 0.0 !!  
    1237           ENDIF 
    1238           IF( med_diag%PDLIMSI3%dgsave   ) THEN 
    1239              CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3  ) 
    1240              pdlimsi3(:,:,: )  = 0.0 !!  
    1241           ENDIF 
    1242  
    1243        ENDIF 
    1244        !! lk_iomput                                    
    1245        !! 
     236      !!------------------------------------------------------------------ 
     237      !! Initialise arrays to zero and set up arrays for diagnostics 
     238      !!------------------------------------------------------------------ 
     239      CALL bio_medusa_init( kt ) 
     240 
    1246241# if defined key_axy_nancheck 
    1247        DO jn = 1,jptra 
     242       DO jn = jp_msa0,jp_msa1 
    1248243          !! fq0 = MINVAL(trn(:,:,:,jn)) 
    1249244          !! fq1 = MAXVAL(trn(:,:,:,jn)) 
    1250245          fq2 = SUM(trn(:,:,:,jn)) 
    1251           !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK', & 
    1252           !! &        kt, jn, fq0, fq1, fq2 
    1253           !! AXY (30/01/14): much to our surprise, the next line doesn't work on HECTOR 
    1254           !!                 and has been replaced here with a specialist routine 
     246          !! if (lwp) write (numout,'(a,2i6,3(1x,1pe15.5))') 'NAN-CHECK',     & 
     247          !!                kt, jn, fq0, fq1, fq2 
     248          !! AXY (30/01/14): much to our surprise, the next line doesn't  
     249          !!                 work on HECTOR and has been replaced here with  
     250          !!                 a specialist routine 
    1255251          !! if (fq2 /= fq2 ) then 
    1256252          if ( ieee_is_nan( fq2 ) ) then 
    1257253             !! there's a NaN here 
    1258              if (lwp) write(numout,*) 'NAN detected in field', jn, 'at time', kt, 'at position:' 
     254             if (lwp) write(numout,*) 'NAN detected in field', jn,           & 
     255                                      'at time', kt, 'at position:' 
    1259256             DO jk = 1,jpk 
    1260257                DO jj = 1,jpj 
    … …  
    1263260                      !! if (trn(ji,jj,jk,jn) /= trn(ji,jj,jk,jn)) then 
    1264261                      if ( ieee_is_nan( trn(ji,jj,jk,jn) ) ) then 
    1265                          if (lwp) write (numout,'(a,1pe12.2,4i6)') 'NAN-CHECK', & 
    1266                          &        tmask(ji,jj,jk), ji, jj, jk, jn 
     262                         if (lwp) write (numout,'(a,1pe12.2,4i6)')           & 
     263                            'NAN-CHECK', tmask(ji,jj,jk), ji, jj, jk, jn 
    1267264                      endif 
    1268265                   enddo 
    … …  
    1276273 
    1277274# if defined key_debug_medusa 
    1278       IF (lwp) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 
     275      IF (lwp) write (numout,*)                                              & 
     276                     'trc_bio_medusa: variables initialised and checked' 
    1279277      CALL flush(numout) 
    1280278# endif  
    1281279 
    1282280# if defined key_roam 
    1283       !!---------------------------------------------------------------------- 
     281      !!------------------------------------------------------------------ 
    1284282      !! calculate atmospheric pCO2 
    1285       !!---------------------------------------------------------------------- 
     283      !!------------------------------------------------------------------ 
    1286284      !! 
    1287285      !! what's atmospheric pCO2 doing? (data start in 1859) 
    … …  
    1290288      if (iyr1 .le. 1) then 
    1291289         !! before 1860 
    1292          f_xco2a = hist_pco2(1) 
     290         f_xco2a(:,:) = hist_pco2(1) 
    1293291      elseif (iyr2 .ge. 242) then 
    1294292         !! after 2099 
    1295          f_xco2a = hist_pco2(242) 
     293         f_xco2a(:,:) = hist_pco2(242) 
    1296294      else 
    1297295         !! just right 
    … …  
    1301299         !! AXY (14/06/12): tweaked to make more sense (and be correct) 
    1302300#  if defined key_bs_axy_yrlen 
    1303          fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  !! bugfix: for 360d year with HadGEM2-ES forcing 
     301         !! bugfix: for 360d year with HadGEM2-ES forcing 
     302         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0   
    1304303#  else 
    1305          fq3 = (real(nday_year) - 1.0 + fq2) / 365.0  !! original use of 365 days (not accounting for leap year or 360d year) 
     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 
    1306307#  endif 
    1307308         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
    1308          f_xco2a = fq4 
     309         f_xco2a(:,:) = fq4 
    1309310      endif 
    1310311#  if defined key_axy_pi_co2 
    1311       !! f_xco2a = 284.725       !! CMIP5 pre-industrial pCO2 
    1312       f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2 
     312      !! OCMIP pre-industrial pCO2 
     313      !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2 
     314      f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2  
    1313315#  endif 
    1314316      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    … …  
    1320322      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2 
    1321323      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3 
    1322       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a 
     324      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1) 
    1323325# endif 
    1324326 
    … …  
    1338340      !!============================= 
    1339341      !! Jpalm -- 07-10-2016 -- need to change carb-chem frequency call : 
    1340       !!          we don't want to call on the first time-step of all run submission,  
    1341       !!          but only on the very first time-step, and then every month 
    1342       !!          So we call on nittrc000 if not restarted run,  
    1343       !!          else if one month after last call. 
    1344       !!          assume one month is 30d --> 3600*24*30 : 2592000s 
    1345       !!          try to call carb-chem at 1st month's tm-stp : x * 30d + 1*rdt(i.e: mod = rdt)    
     342      !!          we don't want to call on the first time-step of all run  
     343      !!          submission, but only on the very first time-step, and  
     344      !!          then every month. So we call on nittrc000 if not  
     345      !!          restarted run, else if one month after last call. 
     346      !!          Assume one month is 30d --> 3600*24*30 : 2592000s 
     347      !!          try to call carb-chem at 1st month's tm-stp :  
     348      !!          x * 30d + 1*rdt(i.e: mod = rdt)    
    1346349      !!          ++ need to pass carb-chem output var through restarts 
    1347       !! We want this to be start of month or if starting afresh from  
    1348       !! climatology - marc 20/6/17 
    1349350      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        & 
    1350351           ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN 
    1351          !!---------------------------------------------------------------------- 
     352         !!--------------------------------------------------------------- 
    1352353         !! Calculate the carbonate chemistry for the whole ocean on the first 
    1353354         !! simulation timestep and every month subsequently; the resulting 3D 
    1354355         !! field of omega calcite is used to determine the depth of the CCD 
    1355          !!---------------------------------------------------------------------- 
    1356          !! 
    1357          IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt 
    1358          CALL flush(numout) 
    1359          !! blank flags 
    1360          i2_omcal(:,:) = 0 
    1361          i2_omarg(:,:) = 0 
    1362          !! loop over 3D space 
    1363          DO jk = 1,jpk 
    1364             DO jj = 2,jpjm1 
    1365                DO ji = 2,jpim1 
    1366                   !! OPEN wet point IF..THEN loop 
    1367                   if (tmask(ji,jj,jk).eq.1) then 
    1368                      IF (lk_oasis) THEN 
    1369                         f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling 
    1370                      ENDIF 
    1371                      !! do carbonate chemistry 
    1372                      !! 
    1373                      fdep2 = fsdept(ji,jj,jk)           !! set up level midpoint 
    1374                      !! AXY (28/11/16): local seafloor depth 
    1375                      !!                 previously mbathy(ji,jj) - 1, now mbathy(ji,jj) 
    1376                      jmbathy = mbathy(ji,jj) 
    1377                      !! 
    1378                      !! set up required state variables 
    1379                      zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
    1380                      zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 
    1381                      ztmp = tsn(ji,jj,jk,jp_tem)        !! temperature 
    1382                      zsal = tsn(ji,jj,jk,jp_sal)        !! salinity 
    1383 #  if defined key_mocsy 
    1384                      zsil = max(0.,trn(ji,jj,jk,jpsil))        !! silicic acid 
    1385                      zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
    1386 #  endif 
    1387            !! 
    1388            !! AXY (28/02/14): check input fields 
    1389            if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then 
    1390                         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', & 
    1391                         tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',    & 
    1392                         ji, ',', jj, ',', jk, ') at time', kt 
    1393          IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', & 
    1394          tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    1395                         ztmp = tsb(ji,jj,jk,jp_tem)     !! temperature 
    1396                      endif 
    1397            if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then 
    1398                         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', & 
    1399                         tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    & 
    1400                         ji, ',', jj, ',', jk, ') at time', kt 
    1401                      endif 
    1402                      !! 
    1403                      !! blank input variables not used at this stage (they relate to air-sea flux) 
    1404                      f_kw660 = 1.0 
    1405                      f_pp0   = 1.0 
    1406                      !! 
    1407                      !! calculate carbonate chemistry at grid cell midpoint 
    1408 #  if defined key_mocsy 
    1409                      !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
    1410                      !!                 chemistry package 
    1411                      CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,         &    ! inputs 
    1412                      f_pp0, fdep2, gphit(ji,jj), f_kw660, f_xco2a, 1,                  &    ! inputs 
    1413                      f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),   &    ! outputs 
    1414                      f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,             &    ! outputs 
    1415                      f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,                &    ! outputs 
    1416                      f_co2starair, f_co2flux, f_dpco2 )                                     ! outputs 
    1417                      !! 
    1418                      f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 
    1419                      f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg 
    1420                      f_dcf  = f_rhosw 
    1421 #  else 
    1422                      !! AXY (22/06/15): use old PML carbonate chemistry package (the 
    1423                      !!                 MEDUSA-2 default) 
    1424                      CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660,      &    ! inputs 
    1425                      f_xco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),   &    ! outputs 
    1426                      f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters)      ! outputs 
    1427                      !!  
    1428                      !! AXY (28/02/14): check output fields 
    1429                      if (iters .eq. 25) then 
    1430                         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', & 
    1431                         iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
    1432                      endif 
    1433 #  endif 
    1434                      !! 
    1435                      !! store 3D outputs 
    1436                      f3_pH(ji,jj,jk)    = f_ph 
    1437                      f3_h2co3(ji,jj,jk) = f_h2co3 
    1438                      f3_hco3(ji,jj,jk)  = f_hco3 
    1439                      f3_co3(ji,jj,jk)   = f_co3 
    1440                      f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 
    1441                      f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 
    1442                      !! 
    1443                      !! CCD calculation: calcite 
    1444                      if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then 
    1445                         if (jk .eq. 1) then 
    1446                            f2_ccd_cal(ji,jj) = fdep2 
    1447                         else 
    1448                            fq0 = f3_omcal(ji,jj,jk-1) - f_omcal(ji,jj) 
    1449                            fq1 = f3_omcal(ji,jj,jk-1) - 1.0 
    1450                            fq2 = fq1 / (fq0 + tiny(fq0)) 
    1451                            fq3 = fdep2 - fsdept(ji,jj,jk-1) 
    1452                            fq4 = fq2 * fq3 
    1453                            f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk-1) + fq4 
    1454                         endif 
    1455                         i2_omcal(ji,jj)   = 1 
    1456                      endif 
    1457                      if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
    1458                         !! reached seafloor and still no dissolution; set to seafloor (W-point) 
    1459                         f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1) 
    1460                         i2_omcal(ji,jj)   = 1 
    1461                      endif 
    1462                      !! 
    1463                      !! CCD calculation: aragonite 
    1464                      if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then 
    1465                         if (jk .eq. 1) then 
    1466                            f2_ccd_arg(ji,jj) = fdep2 
    1467                         else 
    1468                            fq0 = f3_omarg(ji,jj,jk-1) - f_omarg(ji,jj) 
    1469                            fq1 = f3_omarg(ji,jj,jk-1) - 1.0 
    1470                            fq2 = fq1 / (fq0 + tiny(fq0)) 
    1471                            fq3 = fdep2 - fsdept(ji,jj,jk-1) 
    1472                            fq4 = fq2 * fq3 
    1473                            f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk-1) + fq4 
    1474                         endif 
    1475                         i2_omarg(ji,jj)   = 1 
    1476                      endif 
    1477                      if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
    1478                         !! reached seafloor and still no dissolution; set to seafloor (W-point) 
    1479                         f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1) 
    1480                         i2_omarg(ji,jj)   = 1 
    1481                      endif 
    1482                   endif 
    1483                ENDDO 
    1484             ENDDO 
    1485          ENDDO 
     356         !!--------------------------------------------------------------- 
     357         CALL carb_chem( kt ) 
     358 
    1486359      ENDIF 
    1487360# endif 
    … …  
    1492365# endif  
    1493366 
    1494       !!---------------------------------------------------------------------- 
     367      !!------------------------------------------------------------------ 
    1495368      !! MEDUSA has unified equation through the water column 
    1496369      !! (Diff. from LOBSTER which has two sets: bio- and non-bio layers)  
    1497370      !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1           
    1498       !!---------------------------------------------------------------------- 
     371      !!------------------------------------------------------------------ 
    1499372      !! 
    1500373      !! NOTE: the ordering of the loops below differs from that of some other 
    … …  
    1512385         !! OPEN horizontal loops 
    1513386         DO jj = 2,jpjm1 
    1514          DO ji = 2,jpim1 
    1515             !! OPEN wet point IF..THEN loop 
    1516             if (tmask(ji,jj,jk).eq.1) then                
    1517                !!====================================================================== 
    1518                !! SETUP LOCAL GRID CELL 
    1519                !!====================================================================== 
    1520                !! 
    1521                !!--------------------------------------------------------------------- 
    1522                !! Some notes on grid vertical structure 
    1523                !! - fsdepw(ji,jj,jk) is the depth of the upper surface of level jk 
    1524                !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk 
    1525                !! - fse3t(ji,jj,jk)  is the thickness of level jk 
    1526                !!--------------------------------------------------------------------- 
    1527                !! 
    1528                !! AXY (11/12/08): set up level thickness 
    1529                fthk  = fse3t(ji,jj,jk) 
    1530                !! AXY (25/02/10): set up level depth (top of level) 
    1531                fdep  = fsdepw(ji,jj,jk) 
    1532                !! AXY (01/03/10): set up level depth (bottom of level) 
    1533                fdep1 = fdep + fthk 
    1534                !! AXY (28/11/16): local seafloor depth 
    1535                !!                 previously mbathy(ji,jj) - 1, now mbathy(ji,jj) 
    1536                jmbathy = mbathy(ji,jj) 
    1537                !! 
    1538                !! set up model tracers 
    1539                !! negative values of state variables are not allowed to 
    1540                !! contribute to the calculated fluxes 
    1541                zchn = max(0.,trn(ji,jj,jk,jpchn)) !! non-diatom chlorophyll 
    1542                zchd = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll 
    1543                zphn = max(0.,trn(ji,jj,jk,jpphn)) !! non-diatoms 
    1544                zphd = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms 
    1545                zpds = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon 
    1546                !! AXY (28/01/10): probably need to take account of chl/biomass connection 
    1547                if (zchn.eq.0.) zphn = 0. 
    1548                if (zchd.eq.0.) zphd = 0. 
    1549                if (zphn.eq.0.) zchn = 0. 
    1550                if (zphd.eq.0.) zchd = 0. 
    1551           !! AXY (23/01/14): duh - why did I forget diatom silicon? 
    1552           if (zpds.eq.0.) zphd = 0. 
    1553           if (zphd.eq.0.) zpds = 0. 
    1554                zzmi = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton 
    1555                zzme = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton 
    1556                zdet = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen 
    1557                zdin = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen 
    1558                zsil = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid 
    1559                zfer = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron" 
    1560 # if defined key_roam 
    1561                zdtc = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon 
    1562                zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
    1563                zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 
    1564                zoxy = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen 
    1565 #  if defined key_axy_carbchem && defined key_mocsy 
    1566                zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
    1567 #  endif 
    1568                !! 
    1569                !! also need physical parameters for gas exchange calculations 
    1570                ztmp = tsn(ji,jj,jk,jp_tem) 
    1571                zsal = tsn(ji,jj,jk,jp_sal) 
    1572                !! 
    1573           !! AXY (28/02/14): check input fields 
    1574                if (ztmp .lt. -3.0 .or. ztmp .gt. 40.0 ) then 
    1575                   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 
    1576                   tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (',    & 
    1577                   ji, ',', jj, ',', jk, ') at time', kt 
    1578         IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', & 
    1579                   tsn(ji,jj,jk,jp_tem), ' -> ', tsb(ji,jj,jk,jp_tem) 
    1580                   ztmp = tsb(ji,jj,jk,jp_tem) !! temperature 
    1581                endif 
    1582                if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then 
    1583                   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 
    1584                   tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (',    & 
    1585                   ji, ',', jj, ',', jk, ') at time', kt 
    1586                endif 
    1587 # else 
    1588                zdtc = zdet * xthetad              !! implicit detrital carbon 
    1589 # endif 
    1590 # if defined key_debug_medusa 
    1591                if (idf.eq.1) then 
    1592                !! AXY (15/01/10) 
    1593                   if (trn(ji,jj,jk,jpdin).lt.0.) then 
    1594                      IF (lwp) write (numout,*) '------------------------------' 
    1595                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR =', trn(ji,jj,jk,jpdin) 
    1596                      IF (lwp) write (numout,*) 'NEGATIVE DIN ERROR @', ji, jj, jk, kt 
    1597                   endif 
    1598                   if (trn(ji,jj,jk,jpsil).lt.0.) then 
    1599                      IF (lwp) write (numout,*) '------------------------------' 
    1600                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR =', trn(ji,jj,jk,jpsil) 
    1601                      IF (lwp) write (numout,*) 'NEGATIVE SIL ERROR @', ji, jj, jk, kt 
    1602                   endif 
    1603 #  if defined key_roam 
    1604                   if (trn(ji,jj,jk,jpdic).lt.0.) then 
    1605                      IF (lwp) write (numout,*) '------------------------------' 
    1606                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR =', trn(ji,jj,jk,jpdic) 
    1607                      IF (lwp) write (numout,*) 'NEGATIVE DIC ERROR @', ji, jj, jk, kt 
    1608                   endif 
    1609                   if (trn(ji,jj,jk,jpalk).lt.0.) then 
    1610                      IF (lwp) write (numout,*) '------------------------------' 
    1611                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR =', trn(ji,jj,jk,jpalk) 
    1612                      IF (lwp) write (numout,*) 'NEGATIVE ALK ERROR @', ji, jj, jk, kt 
    1613                   endif 
    1614                   if (trn(ji,jj,jk,jpoxy).lt.0.) then 
    1615                      IF (lwp) write (numout,*) '------------------------------' 
    1616                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR =', trn(ji,jj,jk,jpoxy) 
    1617                      IF (lwp) write (numout,*) 'NEGATIVE OXY ERROR @', ji, jj, jk, kt 
    1618                   endif 
    1619 #  endif 
    1620                endif 
    1621 # endif 
    1622 # if defined key_debug_medusa 
    1623                !! report state variable values 
    1624                if (idf.eq.1.AND.idfval.eq.1) then 
    1625                   IF (lwp) write (numout,*) '------------------------------' 
    1626                   IF (lwp) write (numout,*) 'fthk(',jk,') = ', fthk 
    1627                   IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn 
    1628                   IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd 
    1629                   IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds 
    1630                   IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi 
    1631                   IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme 
    1632                   IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet 
    1633                   IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin 
    1634                   IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil 
    1635                   IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer 
    1636 #  if defined key_roam 
    1637                   IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc 
    1638                   IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic 
    1639                   IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk 
    1640                   IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy                   
    1641 #  endif 
    1642                endif 
    1643 # endif 
    1644  
    1645 # if defined key_debug_medusa 
    1646                if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 
    1647                   IF (lwp) write (numout,*) '------------------------------' 
    1648                   IF (lwp) write (numout,*) 'dust      = ', dust(ji,jj) 
    1649                endif 
    1650 # endif 
    1651  
    1652                !! sum tracers for inventory checks 
    1653                IF( lk_iomput ) THEN 
    1654                   IF ( med_diag%INVTN%dgsave )   THEN 
    1655                      ftot_n(ji,jj)  = ftot_n(ji,jj) + & 
    1656                              (fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) ) 
    1657                   ENDIF 
    1658                   IF ( med_diag%INVTSI%dgsave )  THEN 
    1659                      ftot_si(ji,jj) = ftot_si(ji,jj) + &  
    1660                              (fthk * ( zpds + zsil ) ) 
    1661                   ENDIF 
    1662                   IF ( med_diag%INVTFE%dgsave )  THEN 
    1663                      ftot_fe(ji,jj) = ftot_fe(ji,jj) + &  
    1664                              (fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) ) 
    1665                   ENDIF 
    1666 # if defined key_roam 
    1667                   IF ( med_diag%INVTC%dgsave )  THEN 
    1668                      ftot_c(ji,jj)  = ftot_c(ji,jj) + &  
    1669                              (fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + & 
    1670                              (xthetazmi * zzmi) + (xthetazme * zzme) + zdtc +   & 
    1671                              zdic ) ) 
    1672                   ENDIF 
    1673                   IF ( med_diag%INVTALK%dgsave ) THEN 
    1674                      ftot_a(ji,jj)  = ftot_a(ji,jj) + (fthk * ( zalk ) ) 
    1675                   ENDIF 
    1676                   IF ( med_diag%INVTO2%dgsave )  THEN 
    1677                      ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fthk * ( zoxy ) ) 
    1678                   ENDIF 
     387            DO ji = 2,jpim1 
     388               !! OPEN wet point IF..THEN loop 
     389               if (tmask(ji,jj,jk) == 1) then                
     390                  !!====================================================== 
     391                  !! SETUP LOCAL GRID CELL 
     392                  !!====================================================== 
    1679393                  !! 
    1680                   !! AXY (10/11/16): CMIP6 diagnostics 
    1681                   IF ( med_diag%INTDISSIC%dgsave ) THEN 
    1682                      intdissic(ji,jj) = intdissic(ji,jj) + (fthk * zdic) 
    1683                   ENDIF 
    1684                   IF ( med_diag%INTDISSIN%dgsave ) THEN 
    1685                      intdissin(ji,jj) = intdissin(ji,jj) + (fthk * zdin) 
    1686                   ENDIF 
    1687                   IF ( med_diag%INTDISSISI%dgsave ) THEN 
    1688                      intdissisi(ji,jj) = intdissisi(ji,jj) + (fthk * zsil) 
    1689                   ENDIF 
    1690                   IF ( med_diag%INTTALK%dgsave ) THEN 
    1691                      inttalk(ji,jj) = inttalk(ji,jj) + (fthk * zalk) 
    1692                   ENDIF 
    1693                   IF ( med_diag%O2min%dgsave ) THEN 
    1694                      if ( zoxy < o2min(ji,jj) ) then 
    1695                         o2min(ji,jj)  = zoxy 
    1696                         IF ( med_diag%ZO2min%dgsave ) THEN 
    1697                            zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint 
    1698                         ENDIF 
    1699                      endif 
    1700                   ENDIF 
    1701 # endif 
    1702                ENDIF 
    1703  
    1704                CALL flush(numout) 
    1705  
    1706                !!====================================================================== 
    1707                !! LOCAL GRID CELL CALCULATIONS 
    1708                !!====================================================================== 
    1709                !! 
    1710 # if defined key_roam 
    1711                if ( jk .eq. 1 ) then 
    1712                   !!---------------------------------------------------------------------- 
    1713                   !! Air-sea gas exchange 
    1714                   !!---------------------------------------------------------------------- 
     394                  !!------------------------------------------------------ 
     395                  !! Some notes on grid vertical structure 
     396                  !! - fsdepw(ji,jj,jk) is the depth of the upper surface of  
     397                  !!   level jk 
     398                  !! - fsde3w(ji,jj,jk) is *approximately* the midpoint of  
     399                  !!   level jk 
     400                  !! - fse3t(ji,jj,jk)  is the thickness of level jk 
     401                  !!------------------------------------------------------ 
    1715402                  !! 
    1716                   !! AXY (17/07/14): zwind_i and zwind_j do not exist in this 
    1717                   !!                 version of NEMO because it does not include 
    1718                   !!                 the SBC changes that our local version has 
    1719                   !!                 for accessing the HadGEM2 forcing; they  
    1720                   !!                 could be added, but an alternative approach 
    1721                   !!                 is to make use of wndm from oce_trc.F90 
    1722                   !!                 which is wind speed at 10m (which is what 
    1723                   !!                 is required here; this may need to be 
    1724                   !!                 revisited when MEDUSA properly interacts 
    1725                   !!                 with UKESM1 physics 
     403                  !! AXY (01/03/10): set up level depth (bottom of level) 
     404                  fdep1(ji,jj) = fsdepw(ji,jj,jk) + fse3t(ji,jj,jk) 
    1726405                  !! 
    1727                   f_wind  = wndm(ji,jj) 
    1728                   IF (lk_oasis) THEN 
    1729                      f_xco2a = PCO2a_in_cpl(ji,jj)        !! use 2D atm xCO2 from atm coupling 
    1730                   ENDIF 
    1731                   !! 
    1732                   !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
    1733                   !!                 in MEDUSA, the gas transfer velocity used in the carbon 
    1734                   !!                 and oxygen cycles has been harmonised and is calculated 
    1735                   !!                 by the same function here; this harmonisation includes 
    1736                   !!                 changes to the PML carbonate chemistry scheme so that 
    1737                   !!                 it too makes use of the same gas transfer velocity; the 
    1738                   !!                 preferred parameterisation of this is Wanninkhof (2014), 
    1739                   !!                 option 7 
    1740                   !! 
    1741 #   if defined key_debug_medusa 
    1742                      IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
    1743                      CALL flush(numout) 
    1744 #   endif 
    1745                   CALL gas_transfer( f_wind, 1, 7, &  ! inputs 
    1746                                      f_kw660 )        ! outputs 
    1747 #   if defined key_debug_medusa 
    1748                      IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
    1749                      CALL flush(numout) 
    1750 #   endif 
    1751                   !! 
    1752                   !! air pressure (atm); ultimately this will use air pressure at the base 
    1753                   !! of the UKESM1 atmosphere  
    1754                   !!                                      
    1755                   f_pp0   = 1.0 
    1756                   !! 
    1757                   !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp    =', ztmp 
    1758                   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
    1759                   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 
    1760                   !! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind  =', f_wind 
    1761                   !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i    =', fr_i(ji,jj) 
    1762                   !! 
    1763 #  if defined key_axy_carbchem 
    1764 #   if defined key_mocsy 
    1765                   !! 
    1766                   !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY-2 carbonate 
    1767                   !!                 chemistry package; note that depth is set to 
    1768                   !!                 zero in this call 
    1769                   CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho,        &  ! inputs 
    1770                   f_pp0, 0.0, gphit(ji,jj), f_kw660, f_xco2a, 1,                   &  ! inputs 
    1771                   f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj),  &  ! outputs 
    1772                   f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut,            &  ! outputs 
    1773                   f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0,               &  ! outputs 
    1774                   f_co2starair, f_co2flux, f_dpco2 )                                  ! outputs 
    1775                   !! 
    1776                   f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 -> umol / kg 
    1777                   f_TALK = (zalk / f_rhosw) * 1000. !  meq / m3 ->  ueq / kg 
    1778                   f_dcf  = f_rhosw 
    1779 #   else                   
    1780                   iters = 0 
    1781                   !! 
    1782                   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP-2, but not) 
    1783                   CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a,  &  ! inputs 
    1784                   f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj),               &  ! outputs 
    1785                   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters )      ! outputs 
    1786                   !! 
    1787                   !! AXY (09/01/14): removed iteration and NaN checks; these have 
    1788                   !!                 been moved to trc_co2_medusa together with a 
    1789                   !!                 fudge that amends erroneous values (this is 
    1790                   !!                 intended to be a temporary fudge!); the 
    1791                   !!                 output warnings are retained here so that 
    1792                   !!                 failure position can be determined 
    1793                   if (iters .eq. 25) then 
    1794                      IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & 
    1795                      iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
    1796                   endif 
    1797 #   endif 
    1798 #  else 
    1799                   !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
    1800                   !!                 quasi-sensible alternatives 
    1801                   f_ph           = 8.1 
    1802                   f_pco2w        = f_xco2a 
    1803                   f_h2co3        = 0.005 * zdic 
    1804                   f_hco3         = 0.865 * zdic 
    1805                   f_co3          = 0.130 * zdic 
    1806                   f_omcal(ji,jj) = 4. 
    1807                   f_omarg(ji,jj) = 2. 
    1808                   f_co2flux      = 0. 
    1809                   f_TDIC         = zdic 
    1810                   f_TALK         = zalk 
    1811                   f_dcf          = 1.026 
    1812                   f_henry        = 1. 
    1813                   !! AXY (23/06/15): add in some extra MOCSY diagnostics 
    1814                   f_fco2w        = f_xco2a 
    1815                   f_BetaD        = 1. 
    1816                   f_rhosw        = 1.026 
    1817                   f_opres        = 0. 
    1818                   f_insitut      = ztmp 
    1819                   f_pco2atm      = f_xco2a 
    1820                   f_fco2atm      = f_xco2a 
    1821                   f_schmidtco2   = 660. 
    1822                   f_kwco2        = 0. 
    1823                   f_K0           = 0. 
    1824                   f_co2starair   = f_xco2a 
    1825                   f_dpco2        = 0. 
    1826 #  endif 
    1827                   !! 
    1828                   !! mmol/m2/s -> mmol/m3/d; correct for sea-ice; divide through by layer thickness 
    1829                   f_co2flux = (1. - fr_i(ji,jj)) * f_co2flux * 86400. / fthk 
    1830                   !! 
    1831                   !! oxygen (O2); OCMIP-2 code 
    1832                   !! AXY (23/06/15): amend input list for oxygen to account for common gas 
    1833                   !!                 transfer velocity 
    1834                   !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk,  &  ! inputs 
    1835                   !! f_kw660, f_o2flux, f_o2sat )                                                      ! outputs 
    1836                   CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy,  &  ! inputs 
    1837                   f_kwo2, f_o2flux, f_o2sat )                                ! outputs 
    1838                   !! 
    1839                   !! mmol/m2/s -> mol/m3/d; correct for sea-ice; divide through by layer thickness 
    1840                   f_o2flux  = (1. - fr_i(ji,jj)) * f_o2flux * 86400. / fthk 
    1841                   !! 
    1842                   !! Jpalm (08-2014) 
    1843                   !! DMS surface concentration calculation 
    1844                   !! initialy added for UKESM1 model. 
    1845                   !! using MET-OFFICE subroutine. 
    1846                   !! DMS module only needs Chl concentration and MLD 
    1847                   !! to get an aproximate value of DMS concentration. 
    1848                   !! air-sea fluxes are calculated by atmospheric chemitry model 
    1849                   !! from atm and oc-surface concentrations. 
    1850                   !! 
    1851                   !! AXY (13/03/15): this is amended to calculate all of the DMS 
    1852                   !!                 estimates examined during UKESM1 (see comments 
    1853                   !!                 in trcdms_medusa.F90) 
    1854                   !! 
    1855                   !! AXY (25/05/17): amended to additionally pass DIN limitation as well as [DIN]; 
    1856                   !!                 accounts for differences in nutrient half-saturations; changes 
    1857                   !!                 also made in trc_dms_medusa; this permits an additional DMS 
    1858                   !!                 calculation while retaining the existing Anderson one 
    1859                   !! 
    1860                   IF (jdms .eq. 1) THEN 
    1861                      !! 
    1862                      !! calculate weighted half-saturation for DIN uptake 
    1863                      dms_wtkn = ((zphn * xnln) + (zphd * xnld)) / (zphn + zphd) 
    1864                      !! 
    1865                      !! feed in correct inputs 
    1866                      if (jdms_input .eq. 0) then 
    1867                         !! use instantaneous inputs 
    1868                         dms_nlim = zdin / (zdin + dms_wtkn) 
    1869                         !! 
    1870                         CALL trc_dms_medusa( zchn, zchd,                           &  ! inputs 
    1871                         hmld(ji,jj), qsr(ji,jj),                                   &  ! inputs 
    1872                         zdin, dms_nlim,                                            &  ! inputs 
    1873                         dms_andr, dms_simo, dms_aran, dms_hall, dms_andm )            ! outputs 
    1874                      else 
    1875                         !! use diel-average inputs 
    1876                         dms_nlim = zn_dms_din(ji,jj) / (zn_dms_din(ji,jj) + dms_wtkn) 
    1877                         !! 
    1878                         CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), &  ! inputs 
    1879                         zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj),                      &  ! inputs 
    1880                         zn_dms_din(ji,jj), dms_nlim,                               &  ! inputs 
    1881                         dms_andr, dms_simo, dms_aran, dms_hall, dms_andm )            ! outputs 
    1882                      endif 
    1883                      !! 
    1884                      !! assign correct output to variable passed to atmosphere 
    1885                      if     (jdms_model .eq. 1) then 
    1886                         dms_surf = dms_andr 
    1887                      elseif (jdms_model .eq. 2) then 
    1888                         dms_surf = dms_simo 
    1889                      elseif (jdms_model .eq. 3) then 
    1890                         dms_surf = dms_aran 
    1891                      elseif (jdms_model .eq. 4) then 
    1892                         dms_surf = dms_hall 
    1893                      elseif (jdms_model .eq. 5) then 
    1894                         dms_surf = dms_andm 
    1895                      endif 
    1896                      !! 
    1897                      !! 2D diag through iom_use 
    1898                      IF( lk_iomput ) THEN 
    1899                        IF( med_diag%DMS_SURF%dgsave ) THEN 
    1900                          dms_surf2d(ji,jj) = dms_surf 
    1901                        ENDIF 
    1902                        IF( med_diag%DMS_ANDR%dgsave ) THEN 
    1903                          dms_andr2d(ji,jj) = dms_andr 
    1904                        ENDIF 
    1905                        IF( med_diag%DMS_SIMO%dgsave ) THEN 
    1906                          dms_simo2d(ji,jj) = dms_simo 
    1907                        ENDIF 
    1908                        IF( med_diag%DMS_ARAN%dgsave ) THEN 
    1909                          dms_aran2d(ji,jj) = dms_aran 
    1910                        ENDIF 
    1911                        IF( med_diag%DMS_HALL%dgsave ) THEN 
    1912                          dms_hall2d(ji,jj) = dms_hall 
    1913                        ENDIF 
    1914                        IF( med_diag%DMS_ANDM%dgsave ) THEN 
    1915                          dms_andm2d(ji,jj) = dms_andm 
    1916                        ENDIF 
    1917 #   if defined key_debug_medusa 
    1918                        IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' 
    1919                      CALL flush(numout) 
    1920 #   endif  
    1921                      ENDIF 
    1922                      !! End iom 
    1923                   ENDIF 
    1924                   !! End DMS Loop 
    1925                   !! 
    1926                   !! store 2D outputs 
    1927                   !! 
    1928                   !! JPALM -- 17-11-16 -- put fgco2 out of diag request 
    1929                   !!                    is needed for coupling; pass through restart 
    1930                   !! IF( med_diag%FGCO2%dgsave ) THEN 
    1931                      !! convert from  mol/m2/day to kg/m2/s 
    1932                      fgco2(ji,jj) = f_co2flux * fthk * CO2flux_conv  !! mmol-C/m3/d -> kg-CO2/m2/s 
    1933                   !! ENDIF 
    1934                   IF ( lk_iomput ) THEN 
    1935                       IF( med_diag%ATM_PCO2%dgsave ) THEN 
    1936                          f_pco2a2d(ji,jj) = f_pco2atm 
    1937                       ENDIF 
    1938                       IF( med_diag%OCN_PCO2%dgsave ) THEN 
    1939                          f_pco2w2d(ji,jj) = f_pco2w 
    1940                       ENDIF 
    1941                       IF( med_diag%CO2FLUX%dgsave ) THEN 
    1942                          f_co2flux2d(ji,jj) = f_co2flux * fthk           !! mmol/m3/d -> mmol/m2/d 
    1943                       ENDIF 
    1944                       IF( med_diag%TCO2%dgsave ) THEN 
    1945                          f_TDIC2d(ji,jj) = f_TDIC 
    1946                       ENDIF 
    1947                       IF( med_diag%TALK%dgsave ) THEN 
    1948                          f_TALK2d(ji,jj) = f_TALK 
    1949                       ENDIF 
    1950                       IF( med_diag%KW660%dgsave ) THEN 
    1951                          f_kw6602d(ji,jj) = f_kw660 
    1952                       ENDIF 
    1953                       IF( med_diag%ATM_PP0%dgsave ) THEN 
    1954                          f_pp02d(ji,jj) = f_pp0 
    1955                       ENDIF 
    1956                       IF( med_diag%O2FLUX%dgsave ) THEN 
    1957                          f_o2flux2d(ji,jj) = f_o2flux 
    1958                       ENDIF 
    1959                       IF( med_diag%O2SAT%dgsave ) THEN 
    1960                          f_o2sat2d(ji,jj) = f_o2sat 
    1961                       ENDIF 
    1962                       !! AXY (24/11/16): add in extra MOCSY diagnostics 
    1963                       IF( med_diag%ATM_XCO2%dgsave ) THEN 
    1964                          f_xco2a_2d(ji,jj) = f_xco2a 
    1965                       ENDIF 
    1966                       IF( med_diag%OCN_FCO2%dgsave ) THEN 
    1967                          f_fco2w_2d(ji,jj) = f_fco2w 
    1968                       ENDIF 
    1969                       IF( med_diag%ATM_FCO2%dgsave ) THEN 
    1970                          f_fco2a_2d(ji,jj) = f_fco2atm 
    1971                       ENDIF 
    1972                       IF( med_diag%OCN_RHOSW%dgsave ) THEN 
    1973                          f_ocnrhosw_2d(ji,jj) = f_rhosw 
    1974                       ENDIF 
    1975                       IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
    1976                          f_ocnschco2_2d(ji,jj) = f_schmidtco2 
    1977                       ENDIF 
    1978                       IF( med_diag%OCN_KWCO2%dgsave ) THEN 
    1979                          f_ocnkwco2_2d(ji,jj) = f_kwco2 
    1980                       ENDIF 
    1981                       IF( med_diag%OCN_K0%dgsave ) THEN 
    1982                          f_ocnk0_2d(ji,jj) = f_K0 
    1983                       ENDIF 
    1984                       IF( med_diag%CO2STARAIR%dgsave ) THEN 
    1985                          f_co2starair_2d(ji,jj) = f_co2starair 
    1986                       ENDIF 
    1987                       IF( med_diag%OCN_DPCO2%dgsave ) THEN 
    1988                          f_ocndpco2_2d(ji,jj) = f_dpco2 
    1989                       ENDIF 
    1990                   ENDIF 
    1991                   !!  
    1992                endif 
    1993                !! End jk = 1 loop within ROAM key  
    1994  
    1995                !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 
    1996                IF ( med_diag%O2SAT3%dgsave ) THEN 
    1997                   call oxy_sato( ztmp, zsal, f_o2sat3 ) 
    1998                   o2sat3(ji, jj, jk) = f_o2sat3 
    1999                ENDIF 
    2000  
    2001 # endif 
    2002  
    2003                if ( jk .eq. 1 ) then 
    2004                   !!---------------------------------------------------------------------- 
    2005                   !! River inputs 
    2006                   !!---------------------------------------------------------------------- 
    2007                   !! 
    2008                   !! runoff comes in as        kg / m2 / s 
    2009                   !! used and written out as   m3 / m2 / d (= m / d) 
    2010                   !! where                     1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 
    2011                   !! 
    2012                   !! AXY (17/07/14): the compiler doesn't like this line for some reason; 
    2013                   !!                 as MEDUSA doesn't even use runoff for riverine inputs, 
    2014                   !!                 a temporary solution is to switch off runoff entirely 
    2015                   !!                 here; again, this change is one of several that will  
    2016                   !!                 need revisiting once MEDUSA has bedded down in UKESM1; 
    2017                   !!                 particularly so if the land scheme provides information 
    2018                   !!                 concerning nutrient fluxes 
    2019                   !! 
    2020                   !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 
    2021                   f_runoff(ji,jj) = 0.0 
    2022                   !! 
    2023                   !! nutrients are added via rivers to the model in one of two ways: 
    2024                   !!   1. via river concentration; i.e. the average nutrient concentration 
    2025                   !!      of a river water is described by a spatial file, and this is 
    2026                   !!      multiplied by runoff to give a nutrient flux 
    2027                   !!   2. via direct river flux; i.e. the average nutrient flux due to 
    2028                   !!      rivers is described by a spatial file, and this is simply applied 
    2029                   !!      as a direct nutrient flux (i.e. it does not relate or respond to 
    2030                   !!      model runoff) 
    2031                   !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 
    2032                   !! alkalinity are derived from continent-scale DIC estimates (Huang et al.,  
    2033                   !! 2012) and some Arctic river alkalinity estimates (Katya?) 
    2034                   !!  
    2035                   !! as of 19/07/12, riverine nutrients can now be spread vertically across  
    2036                   !! several grid cells rather than just poured into the surface box; this 
    2037                   !! block of code is still executed, however, to set up the total amounts 
    2038                   !! of nutrient entering via rivers 
    2039                   !! 
    2040                   !! nitrogen 
    2041                   if (jriver_n .eq. 1) then 
    2042                      !! river concentration specified; use runoff to calculate input 
    2043                      f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 
    2044                   elseif (jriver_n .eq. 2) then 
    2045                      !! river flux specified; independent of runoff 
    2046                      f_riv_n(ji,jj) = riv_n(ji,jj) 
    2047                   endif 
    2048                   !! 
    2049                   !! silicon 
    2050                   if (jriver_si .eq. 1) then 
    2051                      !! river concentration specified; use runoff to calculate input 
    2052                      f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 
    2053                   elseif (jriver_si .eq. 2) then 
    2054                      !! river flux specified; independent of runoff 
    2055                      f_riv_si(ji,jj) = riv_si(ji,jj) 
    2056                   endif 
    2057                   !! 
    2058                   !! carbon 
    2059                   if (jriver_c .eq. 1) then 
    2060                      !! river concentration specified; use runoff to calculate input 
    2061                      f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 
    2062                   elseif (jriver_c .eq. 2) then 
    2063                      !! river flux specified; independent of runoff 
    2064                      f_riv_c(ji,jj) = riv_c(ji,jj) 
    2065                   endif 
    2066                   !! 
    2067                   !! alkalinity 
    2068                   if (jriver_alk .eq. 1) then 
    2069                      !! river concentration specified; use runoff to calculate input 
    2070                      f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 
    2071                   elseif (jriver_alk .eq. 2) then 
    2072                      !! river flux specified; independent of runoff 
    2073                      f_riv_alk(ji,jj) = riv_alk(ji,jj) 
    2074                   endif 
    2075  
    2076                endif 
    2077  
    2078                !!---------------------------------------------------------------------- 
    2079                !! Chlorophyll calculations 
    2080                !!---------------------------------------------------------------------- 
    2081                !! 
    2082                !! non-diatoms 
    2083           if (zphn.GT.rsmall) then 
    2084                   fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn))) 
    2085                   faln    = xaln * fthetan 
    2086                else 
    2087                   fthetan = 0. 
    2088                   faln    = 0. 
    2089                endif 
    2090                !! 
    2091                !! diatoms 
    2092           if (zphd.GT.rsmall) then 
    2093                   fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd))) 
    2094                   fald    = xald * fthetad 
    2095                else 
    2096                   fthetad = 0. 
    2097                   fald    = 0. 
    2098                endif 
    2099  
    2100 # if defined key_debug_medusa 
    2101                !! report biological calculations 
    2102                if (idf.eq.1.AND.idfval.eq.1) then 
    2103                   IF (lwp) write (numout,*) '------------------------------' 
    2104                   IF (lwp) write (numout,*) 'faln(',jk,') = ', faln 
    2105                   IF (lwp) write (numout,*) 'fald(',jk,') = ', fald 
    2106                endif 
    2107 # endif 
    2108  
    2109                !!---------------------------------------------------------------------- 
    2110                !! Phytoplankton light limitation 
    2111                !!---------------------------------------------------------------------- 
    2112                !! 
    2113                !! It is assumed xpar is the depth-averaged (vertical layer) PAR  
    2114                !! Light limitation (check self-shading) in W/m2 
    2115                !! 
    2116                !! Note that there is no temperature dependence in phytoplankton 
    2117                !! growth rate or any other function.  
    2118                !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid 
    2119                !! NaNs in case of Phy==0.   
    2120                !! 
    2121                !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and non-diat:  
    2122                !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012 
    2123                !! 
    2124                !! AXY (16/07/09) 
    2125                !! temperature for new Eppley style phytoplankton growth 
    2126                loc_T   = tsn(ji,jj,jk,jp_tem) 
    2127                fun_T   = 1.066**(1.0 * loc_T) 
    2128                !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for 
    2129                !phytoplankton 
    2130                !!                 growth; remin. unaffected 
    2131                fun_Q10 = jq10**((loc_T - 0.0) / 10.0) 
    2132                if (jphy.eq.1) then 
    2133                   xvpnT = xvpn * fun_T 
    2134                   xvpdT = xvpd * fun_T 
    2135                elseif (jphy.eq.2) then 
    2136                   xvpnT = xvpn * fun_Q10 
    2137                   xvpdT = xvpd * fun_Q10 
    2138                else 
    2139                   xvpnT = xvpn 
    2140                   xvpdT = xvpd 
    2141                endif 
    2142                !! 
    2143                !! non-diatoms 
    2144                fchn1   = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
    2145                if (fchn1.GT.rsmall) then 
    2146                   fchn    = xvpnT / (sqrt(fchn1) + tiny(fchn1)) 
    2147                else 
    2148                   fchn    = 0. 
    2149                endif 
    2150                fjln    = fchn * faln * xpar(ji,jj,jk) !! non-diatom J term 
    2151                fjlim_pn = fjln / xvpnT 
    2152                !! 
    2153                !! diatoms 
    2154                fchd1   = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
    2155                if (fchd1.GT.rsmall) then 
    2156                   fchd    = xvpdT / (sqrt(fchd1) + tiny(fchd1)) 
    2157                else 
    2158                   fchd    = 0. 
    2159                endif 
    2160                fjld    = fchd * fald * xpar(ji,jj,jk) !! diatom J term 
    2161                fjlim_pd = fjld / xvpdT 
    2162        
    2163 # if defined key_debug_medusa 
    2164                !! report phytoplankton light limitation 
    2165                if (idf.eq.1.AND.idfval.eq.1) then 
    2166                   IF (lwp) write (numout,*) '------------------------------' 
    2167                   IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn 
    2168                   IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd 
    2169                   IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln 
    2170                   IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld 
    2171                endif 
    2172 # endif 
    2173  
    2174                !!---------------------------------------------------------------------- 
    2175                !! Phytoplankton nutrient limitation 
    2176                !!---------------------------------------------------------------------- 
    2177                !! 
    2178                !! non-diatoms (N, Fe) 
    2179                fnln = zdin / (zdin + xnln) !! non-diatom Qn term 
    2180                ffln = zfer / (zfer + xfln) !! non-diatom Qf term 
    2181                !! 
    2182                !! diatoms (N, Si, Fe) 
    2183                fnld = zdin / (zdin + xnld) !! diatom Qn term 
    2184                fsld = zsil / (zsil + xsld) !! diatom Qs term 
    2185                ffld = zfer / (zfer + xfld) !! diatom Qf term 
    2186  
    2187 # if defined key_debug_medusa 
    2188                !! report phytoplankton nutrient limitation 
    2189                if (idf.eq.1.AND.idfval.eq.1) then 
    2190                   IF (lwp) write (numout,*) '------------------------------' 
    2191                   IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln 
    2192                   IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld 
    2193                   IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln 
    2194                   IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld 
    2195                   IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld 
    2196                endif 
    2197 # endif 
    2198  
    2199                !!---------------------------------------------------------------------- 
    2200                !! Primary production (non-diatoms) 
    2201                !! (note: still needs multiplying by phytoplankton concentration) 
    2202                !!---------------------------------------------------------------------- 
    2203                !! 
    2204                if (jliebig .eq. 0) then 
    2205                   !! multiplicative nutrient limitation 
    2206                   fpnlim = fnln * ffln 
    2207                elseif (jliebig .eq. 1) then 
    2208                   !! Liebig Law (= most limiting) nutrient limitation 
    2209                   fpnlim = min(fnln, ffln) 
    2210                endif 
    2211                fprn = fjln * fpnlim 
    2212  
    2213                !!---------------------------------------------------------------------- 
    2214                !! Primary production (diatoms) 
    2215                !! (note: still needs multiplying by phytoplankton concentration) 
    2216                !! 
    2217                !! production here is split between nitrogen production and that of 
    2218                !! silicon; depending upon the "intracellular" ratio of Si:N, model 
    2219                !! diatoms will uptake nitrogen/silicon differentially; this borrows 
    2220                !! from the diatom model of Mongin et al. (2006) 
    2221                !!---------------------------------------------------------------------- 
    2222                !! 
    2223                if (jliebig .eq. 0) then 
    2224                   !! multiplicative nutrient limitation 
    2225                   fpdlim = fnld * ffld 
    2226                elseif (jliebig .eq. 1) then 
    2227                   !! Liebig Law (= most limiting) nutrient limitation 
    2228                   fpdlim = min(fnld, ffld) 
    2229                endif 
    2230                !! 
    2231           if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then 
    2232                   !! "intracellular" elemental ratios 
    2233                   ! fsin  = zpds / (zphd + tiny(zphd)) 
    2234                   ! fnsi  = zphd / (zpds + tiny(zpds)) 
    2235                   fsin = 0.0 
    2236                   IF( zphd .GT. rsmall) fsin  = zpds / zphd 
    2237                   fnsi = 0.0 
    2238                   IF( zpds .GT. rsmall) fnsi  = zphd / zpds 
    2239                   !! AXY (23/02/10): these next variables derive from Mongin et al. (2003) 
    2240                   fsin1 = 3.0 * xsin0 !! = 0.6 
    2241                   fnsi1 = 1.0 / fsin1 !! = 1.667 
    2242                   fnsi2 = 1.0 / xsin0 !! = 5.0 
    2243                   !! 
    2244                   !! conditionalities based on ratios 
    2245                   !! nitrogen (and iron and carbon) 
    2246                   if (fsin.le.xsin0) then 
    2247                      fprd  = 0.0 
    2248                      fsld2 = 0.0 
    2249                   elseif (fsin.lt.fsin1) then 
    2250                      fprd  = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim) 
    2251                      fsld2 = xuif * ((fsin - xsin0) / (fsin + tiny(fsin))) 
    2252                   elseif (fsin.ge.fsin1) then 
    2253                      fprd  = (fjld * fpdlim) 
    2254                      fsld2 = 1.0 
    2255                   endif 
    2256                   !! 
    2257                   !! silicon 
    2258                   if (fsin.lt.fnsi1) then 
    2259                      fprds = (fjld * fsld) 
    2260                   elseif (fsin.lt.fnsi2) then 
    2261                      fprds = xuif * ((fnsi - xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld) 
    2262                   else 
    2263                      fprds = 0.0 
    2264                   endif      
    2265                else 
    2266                   fsin  = 0.0 
    2267                   fnsi  = 0.0 
    2268                   fprd  = 0.0 
    2269                   fsld2 = 0.0 
    2270                   fprds = 0.0 
    2271                endif 
    2272  
    2273 # if defined key_debug_medusa 
    2274                !! report phytoplankton growth (including diatom silicon submodel) 
    2275                if (idf.eq.1.AND.idfval.eq.1) then 
    2276                   IF (lwp) write (numout,*) '------------------------------' 
    2277                   IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin 
    2278                   IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi 
    2279                   IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2 
    2280                   IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn 
    2281                   IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd 
    2282                   IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds 
    2283                endif 
    2284 # endif 
    2285  
    2286                !!---------------------------------------------------------------------- 
    2287                !! Mixed layer primary production 
    2288                !! this block calculates the amount of primary production that occurs 
    2289                !! within the upper mixed layer; this allows the separate diagnosis 
    2290                !! of "sub-surface" primary production; it does assume that short- 
    2291                !! term variability in mixed layer depth doesn't mess with things 
    2292                !! though 
    2293                !!---------------------------------------------------------------------- 
    2294                !! 
    2295                if (fdep1.le.hmld(ji,jj)) then 
    2296                   !! this level is entirely in the mixed layer 
    2297                   fq0 = 1.0 
    2298                elseif (fdep.ge.hmld(ji,jj)) then 
    2299                   !! this level is entirely below the mixed layer 
    2300                   fq0 = 0.0 
    2301                else 
    2302                   !! this level straddles the mixed layer 
    2303                   fq0 = (hmld(ji,jj) - fdep) / fthk 
    2304                endif 
    2305                !! 
    2306                fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0) 
    2307                fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0) 
    2308                 
    2309                !!---------------------------------------------------------------------- 
    2310                !! Vertical Integral -- 
    2311                !!---------------------------------------------------------------------- 
    2312                ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn * fthk)   !! vertical integral non-diatom phytoplankton 
    2313                ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd * fthk)   !! vertical integral diatom phytoplankton 
    2314                ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk)   !! vertical integral microzooplankton 
    2315                ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk)   !! vertical integral mesozooplankton 
    2316                ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk)   !! vertical integral slow detritus, nitrogen 
    2317                ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk)   !! vertical integral slow detritus, carbon 
    2318                 
    2319                !!---------------------------------------------------------------------- 
    2320                !! More chlorophyll calculations 
    2321                !!---------------------------------------------------------------------- 
    2322                !! 
    2323                !! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk))) 
    2324                !! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk))) 
    2325                frn = (xthetam * fchn * fnln * ffln       ) / (fthetan + tiny(fthetan)) 
    2326                !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid  
    2327                !!   limitation, is used in the following line to regulate chlorophyll  
    2328                !!   growth in a manner that is inconsistent with its use in the regulation  
    2329                !!   of biomass growth; the Mongin term term used in growth is more complex 
    2330                !!   than the simple multiplicative function used below 
    2331                !! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad)) 
    2332                !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to 
    2333                !!   regulate chlorophyll growth 
    2334                frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad)) 
    2335  
    2336 # if defined key_debug_medusa 
    2337                !! report chlorophyll calculations 
    2338                if (idf.eq.1.AND.idfval.eq.1) then 
    2339                   IF (lwp) write (numout,*) '------------------------------' 
    2340                   IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan 
    2341                   IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad 
    2342                   IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn 
    2343                   IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd 
    2344                endif 
    2345 # endif 
    2346  
    2347                !!---------------------------------------------------------------------- 
    2348                !! Zooplankton Grazing  
    2349                !! this code supplements the base grazing model with one that 
    2350                !! considers the C:N ratio of grazed food and balances this against 
    2351                !! the requirements of zooplankton growth; this model is derived  
    2352                !! from that of Anderson & Pondaven (2003) 
    2353                !! 
    2354                !! the current version of the code assumes a fixed C:N ratio for 
    2355                !! detritus (in contrast to Anderson & Pondaven, 2003), though the 
    2356                !! full equations are retained for future extension 
    2357                !!---------------------------------------------------------------------- 
    2358                !! 
    2359                !!---------------------------------------------------------------------- 
    2360                !! Microzooplankton first 
    2361                !!---------------------------------------------------------------------- 
    2362                !! 
    2363                fmi1    = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet) 
    2364                fmi     = xgmi * zzmi / fmi1 
    2365                fgmipn  = fmi * xpmipn * zphn * zphn   !! grazing on non-diatoms 
    2366                fgmid   = fmi * xpmid  * zdet * zdet   !! grazing on detrital nitrogen 
    2367 # if defined key_roam 
    2368                fgmidc  = rsmall !acc 
    2369                IF ( zdet .GT. rsmall ) fgmidc  = (zdtc / (zdet + tiny(zdet))) * fgmid  !! grazing on detrital carbon 
    2370 # else 
    2371                !! AXY (26/11/08): implicit detrital carbon change 
    2372                fgmidc  = xthetad * fgmid              !! grazing on detrital carbon 
    2373 # endif 
    2374                !! 
    2375                !! which translates to these incoming N and C fluxes 
    2376                finmi   = (1.0 - xphi) * (fgmipn + fgmid) 
    2377                ficmi   = (1.0 - xphi) * ((xthetapn * fgmipn) + fgmidc) 
    2378                !! 
    2379                !! the ideal food C:N ratio for microzooplankton 
    2380                !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
    2381                fstarmi = (xbetan * xthetazmi) / (xbetac * xkc) 
    2382                !! 
    2383                !! process these to determine proportioning of grazed N and C 
    2384                !! (since there is no explicit consideration of respiration, 
    2385                !! only growth and excretion are calculated here) 
    2386                fmith   = (ficmi / (finmi + tiny(finmi))) 
    2387                if (fmith.ge.fstarmi) then 
    2388                   fmigrow = xbetan * finmi 
    2389                   fmiexcr = 0.0 
    2390                else 
    2391                   fmigrow = (xbetac * xkc * ficmi) / xthetazmi 
    2392                   fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith))) - ((xbetac * xkc) / xthetazmi)) 
    2393                endif 
    2394 # if defined key_roam 
    2395                fmiresp = (xbetac * ficmi) - (xthetazmi * fmigrow) 
    2396 # endif 
    2397  
    2398 # if defined key_debug_medusa 
    2399                !! report microzooplankton grazing 
    2400                if (idf.eq.1.AND.idfval.eq.1) then 
    2401                   IF (lwp) write (numout,*) '------------------------------' 
    2402                   IF (lwp) write (numout,*) 'fmi1(',jk,')    = ', fmi1 
    2403                   IF (lwp) write (numout,*) 'fmi(',jk,')     = ', fmi 
    2404                   IF (lwp) write (numout,*) 'fgmipn(',jk,')  = ', fgmipn 
    2405                   IF (lwp) write (numout,*) 'fgmid(',jk,')   = ', fgmid 
    2406                   IF (lwp) write (numout,*) 'fgmidc(',jk,')  = ', fgmidc 
    2407                   IF (lwp) write (numout,*) 'finmi(',jk,')   = ', finmi 
    2408                   IF (lwp) write (numout,*) 'ficmi(',jk,')   = ', ficmi 
    2409                   IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi 
    2410                   IF (lwp) write (numout,*) 'fmith(',jk,')   = ', fmith 
    2411                   IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow 
    2412                   IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr 
    2413 #  if defined key_roam 
    2414                   IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp 
    2415 #  endif 
    2416                endif 
    2417 # endif 
    2418  
    2419                !!---------------------------------------------------------------------- 
    2420                !! Mesozooplankton second 
    2421                !!---------------------------------------------------------------------- 
    2422                !! 
    2423                fme1    = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + &  
    2424                          (xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet) 
    2425                fme     = xgme * zzme / fme1 
    2426                fgmepn  = fme * xpmepn  * zphn * zphn  !! grazing on non-diatoms 
    2427                fgmepd  = fme * xpmepd  * zphd * zphd  !! grazing on diatoms 
    2428                fgmepds = fsin * fgmepd                !! grazing on diatom silicon 
    2429                fgmezmi = fme * xpmezmi * zzmi * zzmi  !! grazing on microzooplankton 
    2430                fgmed   = fme * xpmed   * zdet * zdet  !! grazing on detrital nitrogen 
    2431 # if defined key_roam 
    2432                fgmedc  = rsmall !acc 
    2433                IF ( zdet .GT. rsmall ) fgmedc  = (zdtc / (zdet + tiny(zdet))) * fgmed  !! grazing on detrital carbon 
    2434 # else 
    2435                !! AXY (26/11/08): implicit detrital carbon change 
    2436                fgmedc  = xthetad * fgmed              !! grazing on detrital carbon 
    2437 # endif 
    2438                !! 
    2439                !! which translates to these incoming N and C fluxes 
    2440                finme   = (1.0 - xphi) * (fgmepn + fgmepd + fgmezmi + fgmed) 
    2441                ficme   = (1.0 - xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & 
    2442                         (xthetazmi * fgmezmi) + fgmedc) 
    2443                !! 
    2444                !! the ideal food C:N ratio for mesozooplankton 
    2445                !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
    2446                fstarme = (xbetan * xthetazme) / (xbetac * xkc) 
    2447                !! 
    2448                !! process these to determine proportioning of grazed N and C 
    2449                !! (since there is no explicit consideration of respiration, 
    2450                !! only growth and excretion are calculated here) 
    2451                fmeth   = (ficme / (finme + tiny(finme))) 
    2452                if (fmeth.ge.fstarme) then 
    2453                   fmegrow = xbetan * finme 
    2454                   fmeexcr = 0.0 
    2455                else 
    2456                   fmegrow = (xbetac * xkc * ficme) / xthetazme 
    2457                   fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth))) - ((xbetac * xkc) / xthetazme)) 
    2458                endif 
    2459 # if defined key_roam 
    2460                fmeresp = (xbetac * ficme) - (xthetazme * fmegrow) 
    2461 # endif 
    2462  
    2463 # if defined key_debug_medusa 
    2464                !! report mesozooplankton grazing 
    2465                if (idf.eq.1.AND.idfval.eq.1) then 
    2466                   IF (lwp) write (numout,*) '------------------------------' 
    2467                   IF (lwp) write (numout,*) 'fme1(',jk,')    = ', fme1 
    2468                   IF (lwp) write (numout,*) 'fme(',jk,')     = ', fme 
    2469                   IF (lwp) write (numout,*) 'fgmepn(',jk,')  = ', fgmepn 
    2470                   IF (lwp) write (numout,*) 'fgmepd(',jk,')  = ', fgmepd 
    2471                   IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds 
    2472                   IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi 
    2473                   IF (lwp) write (numout,*) 'fgmed(',jk,')   = ', fgmed 
    2474                   IF (lwp) write (numout,*) 'fgmedc(',jk,')  = ', fgmedc 
    2475                   IF (lwp) write (numout,*) 'finme(',jk,')   = ', finme 
    2476                   IF (lwp) write (numout,*) 'ficme(',jk,')   = ', ficme 
    2477                   IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme 
    2478                   IF (lwp) write (numout,*) 'fmeth(',jk,')   = ', fmeth 
    2479                   IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow 
    2480                   IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr 
    2481 #  if defined key_roam 
    2482                   IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp 
    2483 #  endif 
    2484                endif 
    2485 # endif 
    2486  
    2487                fzmi_i(ji,jj)  = fzmi_i(ji,jj)  + fthk * (  & 
    2488                   fgmipn + fgmid ) 
    2489                fzmi_o(ji,jj)  = fzmi_o(ji,jj)  + fthk * (  & 
    2490                   fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0 - xbetan) * finmi) ) 
    2491                fzme_i(ji,jj)  = fzme_i(ji,jj)  + fthk * (  & 
    2492                   fgmepn + fgmepd + fgmezmi + fgmed ) 
    2493                fzme_o(ji,jj)  = fzme_o(ji,jj)  + fthk * (  & 
    2494                   fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0 - xbetan) * finme) ) 
    2495  
    2496                !!---------------------------------------------------------------------- 
    2497                !! Plankton metabolic losses 
    2498                !! Linear loss processes assumed to be metabolic in origin 
    2499                !!---------------------------------------------------------------------- 
    2500                !! 
    2501                fdpn2  = xmetapn  * zphn 
    2502                fdpd2  = xmetapd  * zphd 
    2503                fdpds2 = xmetapd  * zpds 
    2504                fdzmi2 = xmetazmi * zzmi 
    2505                fdzme2 = xmetazme * zzme 
    2506  
    2507                !!---------------------------------------------------------------------- 
    2508                !! Plankton mortality losses 
    2509                !! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced  
    2510                !! to improve performance in gyres 
    2511                !!---------------------------------------------------------------------- 
    2512                !! 
    2513                !! non-diatom phytoplankton 
    2514                if (jmpn.eq.1) fdpn = xmpn * zphn               !! linear 
    2515                if (jmpn.eq.2) fdpn = xmpn * zphn * zphn        !! quadratic 
    2516                if (jmpn.eq.3) fdpn = xmpn * zphn * &           !! hyperbolic 
    2517                   (zphn / (xkphn + zphn)) 
    2518                if (jmpn.eq.4) fdpn = xmpn * zphn * &           !! sigmoid 
    2519                   ((zphn * zphn) / (xkphn + (zphn * zphn))) 
    2520                !! 
    2521                !! diatom phytoplankton 
    2522                if (jmpd.eq.1) fdpd = xmpd * zphd               !! linear 
    2523                if (jmpd.eq.2) fdpd = xmpd * zphd * zphd        !! quadratic 
    2524                if (jmpd.eq.3) fdpd = xmpd * zphd * &           !! hyperbolic 
    2525                   (zphd / (xkphd + zphd)) 
    2526                if (jmpd.eq.4) fdpd = xmpd * zphd * &           !! sigmoid 
    2527                   ((zphd * zphd) / (xkphd + (zphd * zphd))) 
    2528                fdpds = fdpd * fsin 
    2529                !! 
    2530                !! microzooplankton 
    2531                if (jmzmi.eq.1) fdzmi = xmzmi * zzmi            !! linear 
    2532                if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi     !! quadratic 
    2533                if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * &        !! hyperbolic 
    2534                   (zzmi / (xkzmi + zzmi)) 
    2535                if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * &        !! sigmoid 
    2536                   ((zzmi * zzmi) / (xkzmi + (zzmi * zzmi))) 
    2537                !! 
    2538                !! mesozooplankton 
    2539                if (jmzme.eq.1) fdzme = xmzme * zzme            !! linear 
    2540                if (jmzme.eq.2) fdzme = xmzme * zzme * zzme     !! quadratic 
    2541                if (jmzme.eq.3) fdzme = xmzme * zzme * &        !! hyperbolic 
    2542                   (zzme / (xkzme + zzme)) 
    2543                if (jmzme.eq.4) fdzme = xmzme * zzme * &        !! sigmoid 
    2544                   ((zzme * zzme) / (xkzme + (zzme * zzme))) 
    2545  
    2546                !!---------------------------------------------------------------------- 
    2547                !! Detritus remineralisation 
    2548                !! Constant or temperature-dependent 
    2549                !!---------------------------------------------------------------------- 
    2550                !! 
    2551                if (jmd.eq.1) then 
    2552                   !! temperature-dependent 
    2553                   fdd  = xmd  * fun_T * zdet 
    2554 # if defined key_roam 
    2555                   fddc = xmdc * fun_T * zdtc 
    2556 # endif 
    2557                elseif (jmd.eq.2) then 
    2558                   !! AXY (16/05/13): add in Q10-based parameterisation (def in nmlst) 
    2559                   !! temperature-dependent 
    2560                   fdd  = xmd  * fun_Q10 * zdet 
    2561 # if defined key_roam 
    2562                   fddc = xmdc * fun_Q10 * zdtc 
    2563 # endif 
    2564                else 
    2565                   !! temperature-independent 
    2566                   fdd  = xmd  * zdet 
    2567 # if defined key_roam 
    2568                   fddc = xmdc * zdtc 
    2569 # endif 
    2570                endif 
    2571                !! 
    2572                !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box 
    2573                if ((jk.eq.jmbathy) .and. jsfd.eq.1) then 
    2574                   fdd  = 1.0  * zdet 
    2575 # if defined key_roam 
    2576                   fddc = 1.0  * zdtc 
    2577 # endif 
    2578                endif 
    2579                 
    2580 # if defined key_debug_medusa 
    2581                !! report plankton mortality and remineralisation 
    2582                if (idf.eq.1.AND.idfval.eq.1) then 
    2583                   IF (lwp) write (numout,*) '------------------------------' 
    2584                   IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2 
    2585                   IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2 
    2586                   IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2 
    2587                   IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2 
    2588                   IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2 
    2589                   IF (lwp) write (numout,*) 'fdpn(',jk,')  = ', fdpn 
    2590                   IF (lwp) write (numout,*) 'fdpd(',jk,')  = ', fdpd 
    2591                   IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds 
    2592                   IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi 
    2593                   IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme 
    2594                   IF (lwp) write (numout,*) 'fdd(',jk,')   = ', fdd 
    2595 #  if defined key_roam 
    2596                   IF (lwp) write (numout,*) 'fddc(',jk,')  = ', fddc 
    2597 #  endif 
    2598                endif 
    2599 # endif 
    2600  
    2601                !!---------------------------------------------------------------------- 
    2602                !! Detritus addition to benthos 
    2603                !! If activated, slow detritus in the bottom box will enter the 
    2604                !! benthic pool 
    2605                !!---------------------------------------------------------------------- 
    2606                !! 
    2607                if ((jk.eq.jmbathy) .and. jorgben.eq.1) then 
    2608                   !! this is the BOTTOM OCEAN BOX -> into the benthic pool! 
    2609                   !! 
    2610                   f_sbenin_n(ji,jj)  = (zdet * vsed * 86400.) 
    2611                   f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn) 
    2612 # if defined key_roam 
    2613                   f_sbenin_c(ji,jj)  = (zdtc * vsed * 86400.) 
    2614 # else 
    2615                   f_sbenin_c(ji,jj)  = (zdet * vsed * 86400. * xthetad) 
    2616 # endif 
    2617                endif 
    2618  
    2619                !!---------------------------------------------------------------------- 
    2620                !! Iron chemistry and fractionation 
    2621                !! following the Parekh et al. (2004) scheme adopted by the Met. 
    2622                !! Office, Medusa models total iron but considers "free" and 
    2623                !! ligand-bound forms for the purposes of scavenging (only "free" 
    2624                !! iron can be scavenged 
    2625                !!---------------------------------------------------------------------- 
    2626                !! 
    2627                !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 
    2628                xFeT        = zfer * 1.e3 
    2629                !! 
    2630                !! calculate fractionation (based on Diat-HadOCC; in turn based on Parekh et al., 2004) 
    2631                xb_coef_tmp = xk_FeL * (xLgT - xFeT) - 1.0 
    2632                xb2M4ac     = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0) 
    2633                !! 
    2634                !! "free" ligand concentration 
    2635                xLgF        = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL 
    2636                !! 
    2637                !! ligand-bound iron concentration 
    2638                xFeL        = xLgT - xLgF 
    2639                !! 
    2640                !! "free" iron concentration (and convert to mmol Fe / m3) 
    2641                xFeF        = (xFeT - xFeL) * 1.e-3 
    2642                xFree(ji,jj)= xFeF / (zfer + tiny(zfer)) 
    2643                !! 
    2644                !! scavenging of iron (multiple schemes); I'm only really happy with the  
    2645                !! first one at the moment - the others involve assumptions (sometimes 
    2646                !! guessed at by me) that are potentially questionable 
    2647                !! 
    2648                if (jiron.eq.1) then 
    2649                   !!---------------------------------------------------------------------- 
    2650                   !! Scheme 1: Dutkiewicz et al. (2005) 
    2651                   !! This scheme includes a single scavenging term based solely on a 
    2652                   !! fixed rate and the availablility of "free" iron 
    2653                   !!---------------------------------------------------------------------- 
    2654                   !! 
    2655                   ffescav     = xk_sc_Fe * xFeF                     ! = mmol/m3/d 
    2656                   !! 
    2657                   !!---------------------------------------------------------------------- 
    2658                   !! 
    2659                   !! Mick's code contains a further (optional) implicit "scavenging" of  
    2660                   !! iron that sets an upper bound on "free" iron concentration, and  
    2661                   !! essentially caps the concentration of total iron as xFeL + "free"  
    2662                   !! iron; since the former is constrained by a fixed total ligand  
    2663                   !! concentration (= 1.0 umol/m3), and the latter isn't allowed above  
    2664                   !! this upper bound, total iron is constrained to a maximum of ... 
    2665                   !! 
    2666                   !!    xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3 
    2667                   !!  
    2668                   !! In Mick's code, the actual value of total iron is reset to this 
    2669                   !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't 
    2670                   !! our favoured approach to tracer updating here (not least because 
    2671                   !! of the leapfrog), so here the amount scavenged is augmented by an 
    2672                   !! additional amount that serves to drag total iron back towards that 
    2673                   !! expected from this limitation on iron concentration ... 
    2674                   !! 
    2675                   xmaxFeF     = min((xFeF * 1.e3), 0.3)             ! = umol/m3 
    2676                   !! 
    2677                   !! Here, the difference between current total Fe and (FeL + Fe') is 
    2678                   !! calculated and added to the scavenging flux already calculated 
    2679                   !! above ... 
    2680                   !! 
    2681                   fdeltaFe    = (xFeT - (xFeL + xmaxFeF)) * 1.e-3   ! = mmol/m3 
    2682                   !! 
    2683                   !! This assumes that the "excess" iron is dissipated with a time- 
    2684                   !! scale of 1 day; seems reasonable to me ... (famous last words) 
    2685                   !! 
    2686                   ffescav     = ffescav + fdeltaFe                  ! = mmol/m3/d 
    2687                   !! 
    2688 # if defined key_deep_fe_fix 
    2689                   !! AXY (17/01/13) 
    2690                   !! stop scavenging for iron concentrations below 0.5 umol / m3 
    2691                   !! at depths greater than 1000 m; this aims to end MEDUSA's 
    2692                   !! continual loss of iron at depth without impacting things 
    2693                   !! at the surface too much; the justification for this is that 
    2694                   !! it appears to be what Mick Follows et al. do in their work 
    2695                   !! (as evidenced by the iron initial condition they supplied 
    2696                   !! me with); to be honest, it looks like Follow et al. do this 
    2697                   !! at shallower depths than 1000 m, but I'll stick with this 
    2698                   !! for now; I suspect that this seemingly arbitrary approach 
    2699                   !! effectively "parameterises" the particle-based scavenging 
    2700                   !! rates that other models use (i.e. at depth there are no 
    2701                   !! sinking particles, so scavenging stops); it might be fun 
    2702                   !! justifying this in a paper though! 
    2703                   !! 
    2704                   if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then 
    2705                      ffescav = 0. 
    2706                   endif 
    2707 # endif 
    2708                   !! 
    2709                elseif (jiron.eq.2) then 
    2710                   !!---------------------------------------------------------------------- 
    2711                   !! Scheme 2: Moore et al. (2004) 
    2712                   !! This scheme includes a single scavenging term that accounts for 
    2713                   !! both suspended and sinking particles in the water column; this 
    2714                   !! term scavenges total iron rather than "free" iron 
    2715                   !!---------------------------------------------------------------------- 
    2716                   !! 
    2717                   !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 
    2718                   xFeT        = zfer * 1.e3 
    2719                   !! 
    2720                   !! this has a base scavenging rate (12% / y) which is modified by local 
    2721                   !! particle concentration and sinking flux (and dust - but I'm ignoring 
    2722                   !! that here for now) and which is accelerated when Fe concentration gets 
    2723                   !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as concentrations 
    2724                   !! below 0.4 nM (= 0.4 umol/m3 = 0.0004 mmol/m3) 
    2725                   !! 
    2726                   !! base scavenging rate (0.12 / y) 
    2727                   fbase_scav = 0.12 / 365.25 
    2728                   !! 
    2729                   !! calculate sinking particle part of scaling factor 
    2730                   !! this takes local fast sinking carbon (mmol C / m2 / d) and 
    2731                   !! gets it into nmol C / cm3 / s ("rdt" below is the number of seconds in 
    2732                   !! a model timestep) 
    2733                   !! 
    2734                   !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.) 
    2735                   !! 
    2736                   !! ... actually, re-reading Moore et al.'s equations, it looks like he uses 
    2737                   !! his sinking flux directly, without scaling it by time-step or anything, 
    2738                   !! so I'll copy this here ... 
    2739                   !! 
    2740                   fscal_sink = ffastc(ji,jj) * 1.e2 
    2741                   !! 
    2742                   !! calculate particle part of scaling factor 
    2743                   !! this totals up the carbon in suspended particles (Pn, Pd, Zmi, Zme, D), 
    2744                   !! which comes out in mmol C / m3 (= nmol C / cm3), and then multiplies it 
    2745                   !! by a magic factor, 0.002, to get it into nmol C / cm2 / s 
    2746                   !! 
    2747                   fscal_part = ((xthetapn * zphn) + (xthetapd * zphd) + (xthetazmi * zzmi) + & 
    2748                   (xthetazme * zzme) + (xthetad * zdet)) * 0.002 
    2749                   !! 
    2750                   !! calculate scaling factor for base scavenging rate 
    2751                   !! this uses the (now correctly scaled) sinking flux and standing 
    2752                   !! particle concentration, divides through by some sort of reference 
    2753                   !! value (= 0.0066 nmol C / cm2 / s) and then uses this, or not if its 
    2754                   !! too high, to rescale the base scavenging rate 
    2755                   !! 
    2756                   fscal_scav = fbase_scav * min(((fscal_sink + fscal_part) / 0.0066), 4.0) 
    2757                   !! 
    2758                   !! the resulting scavenging rate is then scaled further according to the 
    2759                   !! local iron concentration (i.e. diminished in low iron regions; enhanced 
    2760                   !! in high iron regions; less alone in intermediate iron regions) 
    2761                   !! 
    2762                   if (xFeT.lt.0.4) then 
    2763                      !! 
    2764                      !! low iron region 
    2765                      !! 
    2766                      fscal_scav = fscal_scav * (xFeT / 0.4) 
    2767                      !! 
    2768                   elseif (xFeT.gt.0.6) then 
    2769                      !! 
    2770                      !! high iron region 
    2771                      !! 
    2772                      fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4)) 
    2773                      !! 
    2774                   else 
    2775                      !! 
    2776                      !! intermediate iron region: do nothing 
    2777                      !! 
    2778                   endif 
    2779                   !!  
    2780                   !! apply the calculated scavenging rate ... 
    2781                   !! 
    2782                   ffescav = fscal_scav * zfer 
    2783                   !! 
    2784                elseif (jiron.eq.3) then 
    2785                   !!---------------------------------------------------------------------- 
    2786                   !! Scheme 3: Moore et al. (2008) 
    2787                   !! This scheme includes a single scavenging term that accounts for 
    2788                   !! sinking particles in the water column, and includes organic C, 
    2789                   !! biogenic opal, calcium carbonate and dust in this (though the 
    2790                   !! latter is ignored here until I work out what units the incoming 
    2791                   !! "dust" flux is in); this term scavenges total iron rather than  
    2792                   !! "free" iron 
    2793                   !!---------------------------------------------------------------------- 
    2794                   !! 
    2795                   !! total iron concentration (mmol Fe / m3 -> umol Fe / m3) 
    2796                   xFeT        = zfer * 1.e3 
    2797                   !! 
    2798                   !! this has a base scavenging rate which is modified by local 
    2799                   !! particle sinking flux (including dust - but I'm ignoring that  
    2800                   !! here for now) and which is accelerated when Fe concentration 
    2801                   !! is > 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as  
    2802                   !! concentrations < 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3) 
    2803                   !! 
    2804                   !! base scavenging rate (Fe_b in paper; units may be wrong there) 
    2805                   fbase_scav = 0.00384 ! (ng)^-1 cm 
    2806                   !! 
    2807                   !! calculate sinking particle part of scaling factor; this converts 
    2808                   !! mmol / m2 / d fluxes of organic carbon, silicon and calcium 
    2809                   !! carbonate into ng / cm2 / s fluxes; it is assumed here that the 
    2810                   !! mass conversions simply consider the mass of the main element 
    2811                   !! (C, Si and Ca) and *not* the mass of the molecules that they are 
    2812                   !! part of; Moore et al. (2008) is unclear on the conversion that 
    2813                   !! should be used 
    2814                   !! 
    2815                   !! milli -> nano; mol -> gram; /m2 -> /cm2; /d -> /s 
    2816                   fscal_csink  = (ffastc(ji,jj)  * 1.e6 * xmassc  * 1.e-4 / 86400.)      ! ng C  / cm2 / s 
    2817                   fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi * 1.e-4 / 86400.)      ! ng Si / cm2 / s 
    2818                   fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca * 1.e-4 / 86400.)      ! ng Ca / cm2 / s 
    2819                   !!  
    2820                   !! sum up these sinking fluxes and convert to ng / cm by dividing 
    2821                   !! through by a sinking rate of 100 m / d = 1.157 cm / s 
    2822                   fscal_sink   = ((fscal_csink * 6.) + fscal_sisink + fscal_casink) / & 
    2823                   (100. * 1.e3 / 86400)                                                  ! ng / cm 
    2824                   !! 
    2825                   !! now calculate the scavenging rate based upon the base rate and 
    2826                   !! this particle flux scaling; according to the published units, 
    2827                   !! the result actually has *no* units, but as it must be expressed 
    2828                   !! per unit time for it to make any sense, I'm assuming a missing 
    2829                   !! "per second" 
    2830                   fscal_scav = fbase_scav * fscal_sink                                   ! / s 
    2831                   !! 
    2832                   !! the resulting scavenging rate is then scaled further according to the 
    2833                   !! local iron concentration (i.e. diminished in low iron regions; enhanced 
    2834                   !! in high iron regions; less alone in intermediate iron regions) 
    2835                   !! 
    2836                   if (xFeT.lt.0.5) then 
    2837                      !! 
    2838                      !! low iron region (0.5 instead of the 0.4 in Moore et al., 2004) 
    2839                      !! 
    2840                      fscal_scav = fscal_scav * (xFeT / 0.5) 
    2841                      !! 
    2842                   elseif (xFeT.gt.0.6) then 
    2843                      !! 
    2844                      !! high iron region (functional form different in Moore et al., 2004) 
    2845                      !! 
    2846                      fscal_scav = fscal_scav + ((xFeT - 0.6) * 0.00904) 
    2847                      !! 
    2848                   else 
    2849                      !! 
    2850                      !! intermediate iron region: do nothing 
    2851                      !! 
    2852                   endif 
    2853                   !!  
    2854                   !! apply the calculated scavenging rate ... 
    2855                   !! 
    2856                   ffescav = fscal_scav * zfer 
    2857                   !! 
    2858                elseif (jiron.eq.4) then 
    2859                   !!---------------------------------------------------------------------- 
    2860                   !! Scheme 4: Galbraith et al. (2010) 
    2861                   !! This scheme includes two scavenging terms, one for organic, 
    2862                   !! particle-based scavenging, and another for inorganic scavenging; 
    2863                   !! both terms scavenge "free" iron only 
    2864                   !!---------------------------------------------------------------------- 
    2865                   !! 
    2866                   !! Galbraith et al. (2010) present a more straightforward outline of  
    2867                   !! the scheme in Parekh et al. (2005) ... 
    2868                   !!  
    2869                   !! sinking particulate carbon available for scavenging 
    2870                   !! this assumes a sinking rate of 100 m / d (Moore & Braucher, 2008), 
    2871                   xCscav1     = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3 
    2872                   !!  
    2873                   !! scale by Honeyman et al. (1981) exponent coefficient 
    2874                   !! multiply by 1.e-3 to express C flux in g C rather than mg C 
    2875                   xCscav2     = (xCscav1 * 1.e-3)**0.58 
    2876                   !! 
    2877                   !! multiply by Galbraith et al. (2010) scavenging rate 
    2878                   xk_org      = 0.5 ! ((g C m/3)^-1) / d 
    2879                   xORGscav    = xk_org * xCscav2 * xFeF 
    2880                   !! 
    2881                   !! Galbraith et al. (2010) also include an inorganic bit ... 
    2882                   !!  
    2883                   !! this occurs at a fixed rate, again based on the availability of 
    2884                   !! "free" iron 
    2885                   !! 
    2886                   !! k_inorg  = 1000 d**-1 nmol Fe**-0.5 kg**-0.5 
    2887                   !! 
    2888                   !! to implement this here, scale xFeF by 1026 to put in units of 
    2889                   !! umol/m3 which approximately equal nmol/kg 
    2890                   !! 
    2891                   xk_inorg    = 1000. ! ((nmol Fe / kg)^1.5) 
    2892                   xINORGscav  = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e-3 
    2893                   !! 
    2894                   !! sum these two terms together 
    2895                   ffescav     = xORGscav + xINORGscav 
    2896                else 
    2897                   !!---------------------------------------------------------------------- 
    2898                   !! No Scheme: you coward! 
    2899                   !! This scheme puts its head in the sand and eskews any decision about 
    2900                   !! how iron is removed from the ocean; prepare to get deluged in iron 
    2901                   !! you fool! 
    2902                   !!---------------------------------------------------------------------- 
    2903                   ffescav     = 0. 
    2904                endif 
    2905  
    2906                !!---------------------------------------------------------------------- 
    2907                !! Other iron cycle processes 
    2908                !!---------------------------------------------------------------------- 
    2909                !! 
    2910                !! aeolian iron deposition 
    2911                if (jk.eq.1) then 
    2912                   !! zirondep   is in mmol-Fe / m2 / day 
    2913                   !! ffetop     is in mmol-dissolved-Fe / m3 / day 
    2914                   ffetop  = zirondep(ji,jj) * xfe_sol / fthk  
    2915                else 
    2916                   ffetop  = 0.0 
    2917                endif 
    2918                !! 
    2919                !! seafloor iron addition 
    2920                !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down 
    2921                !! if (jk.eq.(mbathy(ji,jj)-1).AND.jk.lt.i1100) then 
    2922                if ((jk.eq.jmbathy).AND.jk.le.i0500) then 
    2923                   !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a 
    2924                   !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value 
    2925                   !! but apply it everywhere 
    2926                   !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 1-37) 
    2927                   ffebot  = (xfe_sed / fthk) 
    2928                else 
    2929                   ffebot  = 0.0 
    2930                endif 
    2931  
    2932                !! AXY (16/12/09): remove iron addition/removal processes 
    2933                !! For the purposes of the quarter degree run, the iron cycle is being pegged to the 
    2934                !! initial condition supplied by Mick Follows via restoration with a 30 day period; 
    2935                !! iron addition at the seafloor is still permitted with the idea that this extra 
    2936                !! iron will be removed by the restoration away from the source 
    2937                !! ffescav = 0.0 
    2938                !! ffetop  = 0.0 
    2939                !! ffebot  = 0.0 
    2940  
    2941 # if defined key_debug_medusa 
    2942                !! report miscellaneous calculations 
    2943                if (idf.eq.1.AND.idfval.eq.1) then 
    2944                   IF (lwp) write (numout,*) '------------------------------' 
    2945                   IF (lwp) write (numout,*) 'xfe_sol  = ', xfe_sol 
    2946                   IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass 
    2947                   IF (lwp) write (numout,*) 'ffetop(',jk,')  = ', ffetop 
    2948                   IF (lwp) write (numout,*) 'ffebot(',jk,')  = ', ffebot 
    2949                   IF (lwp) write (numout,*) 'xFree(',jk,')   = ', xFree(ji,jj) 
    2950                   IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav 
    2951                endif 
    2952 # endif 
    2953  
    2954                !!---------------------------------------------------------------------- 
    2955                !! Miscellaneous 
    2956                !!---------------------------------------------------------------------- 
    2957                !! 
    2958                !! diatom frustule dissolution 
    2959                fsdiss  = xsdiss * zpds 
    2960  
    2961 # if defined key_debug_medusa 
    2962                !! report miscellaneous calculations 
    2963                if (idf.eq.1.AND.idfval.eq.1) then 
    2964                   IF (lwp) write (numout,*) '------------------------------' 
    2965                   IF (lwp) write (numout,*) 'fsdiss(',jk,')  = ', fsdiss 
    2966                endif 
    2967 # endif 
    2968  
    2969                !!---------------------------------------------------------------------- 
    2970                !! Slow detritus creation 
    2971                !!---------------------------------------------------------------------- 
    2972                !! this variable integrates the creation of slow sinking detritus 
    2973                !! to allow the split between fast and slow detritus to be  
    2974                !! diagnosed 
    2975                fslown  = fdpn + fdzmi + ((1.0 - xfdfrac1) * fdpd) + & 
    2976                ((1.0 - xfdfrac2) * fdzme) + ((1.0 - xbetan) * (finmi + finme)) 
    2977                !! 
    2978                !! this variable records the slow detrital sinking flux at this 
    2979                !! particular depth; it is used in the output of this flux at 
    2980                !! standard depths in the diagnostic outputs; needs to be 
    2981                !! adjusted from per second to per day because of parameter vsed 
    2982                fslownflux(ji,jj) = zdet * vsed * 86400. 
    2983 # if defined key_roam 
    2984                !! 
    2985                !! and the same for detrital carbon 
    2986                fslowc  = (xthetapn * fdpn) + (xthetazmi * fdzmi) + & 
    2987                (xthetapd * (1.0 - xfdfrac1) * fdpd) + & 
    2988                (xthetazme * (1.0 - xfdfrac2) * fdzme) + & 
    2989                ((1.0 - xbetac) * (ficmi + ficme)) 
    2990                !! 
    2991                !! this variable records the slow detrital sinking flux at this 
    2992                !! particular depth; it is used in the output of this flux at 
    2993                !! standard depths in the diagnostic outputs; needs to be 
    2994                !! adjusted from per second to per day because of parameter vsed 
    2995                fslowcflux(ji,jj) = zdtc * vsed * 86400. 
    2996 # endif 
    2997  
    2998                !!---------------------------------------------------------------------- 
    2999                !! Nutrient regeneration 
    3000                !! this variable integrates total nitrogen regeneration down the 
    3001                !! watercolumn; its value is stored and output as a 2D diagnostic; 
    3002                !! the corresponding dissolution flux of silicon (from sources 
    3003                !! other than fast detritus) is also integrated; note that, 
    3004                !! confusingly, the linear loss terms from plankton compartments 
    3005                !! are labelled as fdX2 when one might have expected fdX or fdX1 
    3006                !!---------------------------------------------------------------------- 
    3007                !! 
    3008                !! nitrogen 
    3009                fregen   = (( (xphi * (fgmipn + fgmid)) +                &  ! messy feeding 
    3010                (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) +           &  ! messy feeding 
    3011                fmiexcr + fmeexcr + fdd +                                &  ! excretion + D remin. 
    3012                fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk)                    ! linear mortality 
    3013                !! 
    3014                !! silicon 
    3015                fregensi = (( fsdiss + ((1.0 - xfdfrac1) * fdpds) +      &  ! dissolution + non-lin. mortality 
    3016                ((1.0 - xfdfrac3) * fgmepds) +                           &  ! egestion by zooplankton 
    3017                fdpds2) * fthk)                                             ! linear mortality 
    3018 # if defined key_roam 
    3019                !! 
    3020                !! carbon 
    3021                fregenc  = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) +  &  ! messy feeding 
    3022                (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) +     &  ! messy feeding 
    3023                (xthetazmi * fgmezmi) + fgmedc)) +                       &  ! messy feeding 
    3024                fmiresp + fmeresp + fddc +                               &  ! respiration + D remin. 
    3025                (xthetapn * fdpn2) + (xthetapd * fdpd2) +                &  ! linear mortality 
    3026                (xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk)        ! linear mortality 
    3027 # endif 
    3028  
    3029                !!---------------------------------------------------------------------- 
    3030                !! Fast-sinking detritus terms 
    3031                !! "local" variables declared so that conservation can be checked; 
    3032                !! the calculated terms are added to the fast-sinking flux later on 
    3033                !! only after the flux entering this level has experienced some 
    3034                !! remineralisation 
    3035                !! note: these fluxes need to be scaled by the level thickness 
    3036                !!---------------------------------------------------------------------- 
    3037                !! 
    3038                !! nitrogen:   diatom and mesozooplankton mortality 
    3039                ftempn         = b0 * ((xfdfrac1 * fdpd)  + (xfdfrac2 * fdzme)) 
    3040                !! 
    3041                !! silicon:    diatom mortality and grazed diatoms 
    3042                ftempsi        = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds)) 
    3043                !! 
    3044                !! iron:       diatom and mesozooplankton mortality 
    3045                ftempfe        = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn) 
    3046                !! 
    3047                !! carbon:     diatom and mesozooplankton mortality 
    3048                ftempc         = b0 * ((xfdfrac1 * xthetapd * fdpd) + &  
    3049                                 (xfdfrac2 * xthetazme * fdzme)) 
    3050                !! 
    3051 # if defined key_roam 
    3052                if (jrratio.eq.0) then 
    3053                   !! CaCO3:      latitudinally-based fraction of total primary production 
    3054                   !!               absolute latitude of current grid cell 
    3055                   flat           = abs(gphit(ji,jj)) 
    3056                   !!               0.10 at equator; 0.02 at pole 
    3057                   fcaco3         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0)) 
    3058                elseif (jrratio.eq.1) then 
    3059                   !! CaCO3:      Ridgwell et al. (2007) submodel, version 1 
    3060                   !!             this uses SURFACE omega calcite to regulate rain ratio 
    3061                   if (f_omcal(ji,jj).ge.1.0) then 
    3062                      fq1 = (f_omcal(ji,jj) - 1.0)**0.81 
    3063                   else 
    3064                      fq1 = 0. 
    3065                   endif 
    3066                   fcaco3 = xridg_r0 * fq1 
    3067                elseif (jrratio.eq.2) then 
    3068                   !! CaCO3:      Ridgwell et al. (2007) submodel, version 2 
    3069                   !!             this uses FULL 3D omega calcite to regulate rain ratio 
    3070                   if (f3_omcal(ji,jj,jk).ge.1.0) then 
    3071                      fq1 = (f3_omcal(ji,jj,jk) - 1.0)**0.81 
    3072                   else 
    3073                      fq1 = 0. 
    3074                   endif 
    3075                   fcaco3 = xridg_r0 * fq1 
    3076                endif 
    3077 # else 
    3078                !! CaCO3:      latitudinally-based fraction of total primary production 
    3079                !!               absolute latitude of current grid cell 
    3080                flat           = abs(gphit(ji,jj)) 
    3081                !!               0.10 at equator; 0.02 at pole 
    3082                fcaco3         = xcaco3a + ((xcaco3b - xcaco3a) * ((90.0 - flat) / 90.0)) 
    3083 # endif 
    3084                !! AXY (09/03/09): convert CaCO3 production from function of  
    3085                !! primary production into a function of fast-sinking material;  
    3086                !! technically, this is what Dunne et al. (2007) do anyway; they  
    3087                !! convert total primary production estimated from surface  
    3088                !! chlorophyll to an export flux for which they apply conversion  
    3089                !! factors to estimate the various elemental fractions (Si, Ca) 
    3090                ftempca        = ftempc * fcaco3 
    3091  
    3092 # if defined key_debug_medusa 
    3093                !! integrate total fast detritus production 
    3094                if (idf.eq.1) then 
    3095                   fifd_n(ji,jj)  = fifd_n(ji,jj)  + (ftempn  * fthk) 
    3096                   fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi * fthk) 
    3097                   fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe * fthk) 
    3098 #  if defined key_roam 
    3099                   fifd_c(ji,jj)  = fifd_c(ji,jj)  + (ftempc  * fthk) 
    3100 #  endif 
    3101                endif 
    3102  
    3103                !! report quantities of fast-sinking detritus for each component 
    3104                if (idf.eq.1.AND.idfval.eq.1) then 
    3105                   IF (lwp) write (numout,*) '------------------------------' 
    3106                   IF (lwp) write (numout,*) 'fdpd(',jk,')    = ', fdpd 
    3107                   IF (lwp) write (numout,*) 'fdzme(',jk,')   = ', fdzme 
    3108                   IF (lwp) write (numout,*) 'ftempn(',jk,')  = ', ftempn 
    3109                   IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi 
    3110                   IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe 
    3111                   IF (lwp) write (numout,*) 'ftempc(',jk,')  = ', ftempc 
    3112                   IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca 
    3113                   IF (lwp) write (numout,*) 'flat(',jk,')    = ', flat 
    3114                   IF (lwp) write (numout,*) 'fcaco3(',jk,')  = ', fcaco3 
    3115                endif 
    3116 # endif 
    3117  
    3118                !!---------------------------------------------------------------------- 
    3119                !! This version of MEDUSA offers a choice of three methods for 
    3120                !! handling the remineralisation of fast detritus.  All three 
    3121                !! do so in broadly the same way: 
    3122                !! 
    3123                !!   1.  Fast detritus is stored as a 2D array                   [ ffastX  ] 
    3124                !!   2.  Fast detritus is added level-by-level                   [ ftempX  ] 
    3125                !!   3.  Fast detritus is not remineralised in the top box       [ freminX ] 
    3126                !!   4.  Remaining fast detritus is remineralised in the bottom  [ fsedX   ] 
    3127                !!       box 
    3128                !! 
    3129                !! The three remineralisation methods are: 
    3130                !!    
    3131                !!   1.  Ballast model (i.e. that published in Yool et al., 2011) 
    3132                !!  (1b. Ballast-sans-ballast model) 
    3133                !!   2.  Martin et al. (1987) 
    3134                !!   3.  Henson et al. (2011) 
    3135                !!  
    3136                !! The first of these couples C, N and Fe remineralisation to 
    3137                !! the remineralisation of particulate Si and CaCO3, but the  
    3138                !! latter two treat remineralisation of C, N, Fe, Si and CaCO3 
    3139                !! completely separately.  At present a switch within the code 
    3140                !! regulates which submodel is used, but this should be moved 
    3141                !! to the namelist file. 
    3142                !!  
    3143                !! The ballast-sans-ballast submodel is an original development 
    3144                !! feature of MEDUSA in which the ballast submodel's general 
    3145                !! framework and parameterisation is used, but in which there 
    3146                !! is no protection of organic material afforded by ballasting 
    3147                !! minerals.  While similar, it is not the same as the Martin  
    3148                !! et al. (1987) submodel. 
    3149                !! 
    3150                !! Since the three submodels behave the same in terms of 
    3151                !! accumulating sinking material and remineralising it all at 
    3152                !! the seafloor, these portions of the code below are common to 
    3153                !! all three. 
    3154                !!---------------------------------------------------------------------- 
    3155  
    3156                if (jexport.eq.1) then 
    3157                   !!====================================================================== 
    3158                   !! BALLAST SUBMODEL 
    3159                   !!====================================================================== 
    3160                   !!  
    3161                   !!---------------------------------------------------------------------- 
    3162                   !! Fast-sinking detritus fluxes, pt. 1: REMINERALISATION 
    3163                   !! aside from explicitly modelled, slow-sinking detritus, the 
    3164                   !! model includes an implicit representation of detrital 
    3165                   !! particles that sink too quickly to be modelled with 
    3166                   !! explicit state variables; this sinking flux is instead 
    3167                   !! instantaneously remineralised down the water column using 
    3168                   !! the version of Armstrong et al. (2002)'s ballast model 
    3169                   !! used by Dunne et al. (2007); the version of this model 
    3170                   !! here considers silicon and calcium carbonate ballast 
    3171                   !! minerals; this section of the code redistributes the fast 
    3172                   !! sinking material generated locally down the water column; 
    3173                   !! this differs from Dunne et al. (2007) in that fast sinking 
    3174                   !! material is distributed at *every* level below that it is 
    3175                   !! generated, rather than at every level below some fixed 
    3176                   !! depth; this scheme is also different in that sinking material  
    3177                   !! generated in one level is aggregated with that generated by 
    3178                   !! shallower levels; this should make the ballast model more 
    3179                   !! self-consistent (famous last words) 
    3180                   !!---------------------------------------------------------------------- 
    3181                   !! 
    3182                   if (jk.eq.1) then 
    3183                      !! this is the SURFACE OCEAN BOX (no remineralisation) 
    3184                      !! 
    3185                      freminc  = 0.0 
    3186                      freminn  = 0.0 
    3187                      freminfe = 0.0 
    3188                      freminsi = 0.0 
    3189                      freminca = 0.0 
    3190                   elseif (jk.le.jmbathy) then 
    3191                      !! this is an OCEAN BOX (remineralise some material) 
    3192                      !! 
    3193                      !! set up CCD depth to be used depending on user choice 
    3194                      if (jocalccd.eq.0) then 
    3195                         !! use default CCD field 
    3196                         fccd_dep = ocal_ccd(ji,jj) 
    3197                      elseif (jocalccd.eq.1) then 
    3198                         !! use calculated CCD field 
    3199                         fccd_dep = f2_ccd_cal(ji,jj) 
    3200                      endif 
    3201                      !! 
    3202                      !! === organic carbon === 
    3203                      fq0      = ffastc(ji,jj)                        !! how much organic C enters this box        (mol) 
    3204                      if (iball.eq.1) then 
    3205                         fq1      = (fq0 * xmassc)                    !! how much it weighs                        (mass) 
    3206                         fq2      = (ffastca(ji,jj) * xmassca)        !! how much CaCO3 enters this box            (mass) 
    3207                         fq3      = (ffastsi(ji,jj) * xmasssi)        !! how much  opal enters this box            (mass) 
    3208                         fq4      = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C                 (mass) 
    3209                         !! this next term is calculated for C but used for N and Fe as well 
    3210                         !! it needs to be protected in case ALL C is protected 
    3211                         if (fq4.lt.fq1) then 
    3212                            fprotf   = (fq4 / (fq1 + tiny(fq1)))      !! protected fraction of total organic C     (non-dim) 
    3213                         else 
    3214                            fprotf   = 1.0                            !! all organic C is protected                (non-dim) 
    3215                         endif 
    3216                         fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic C   (non-dim) 
    3217                         fq6      = (fq0 * fq5)                       !! how much organic C is unprotected         (mol) 
    3218                         fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected C leaves this box    (mol) 
    3219                         fq8      = (fq7 + (fq0 * fprotf))            !! how much total C leaves this box          (mol) 
    3220                         freminc  = (fq0 - fq8) / fthk                !! C remineralisation in this box            (mol) 
    3221                         ffastc(ji,jj) = fq8                           
    3222 # if defined key_debug_medusa 
    3223                         !! report in/out/remin fluxes of carbon for this level 
    3224                            if (idf.eq.1.AND.idfval.eq.1) then 
    3225                               IF (lwp) write (numout,*) '------------------------------' 
    3226                               IF (lwp) write (numout,*) 'totalC(',jk,')  = ', fq1 
    3227                               IF (lwp) write (numout,*) 'prtctC(',jk,')  = ', fq4 
    3228                               IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf 
    3229                               IF (lwp) write (numout,*) '------------------------------' 
    3230                               IF (lwp) write (numout,*) 'IN   C(',jk,')  = ', fq0 
    3231                               IF (lwp) write (numout,*) 'LOST C(',jk,')  = ', freminc * fthk 
    3232                               IF (lwp) write (numout,*) 'OUT  C(',jk,')  = ', fq8 
    3233                               IF (lwp) write (numout,*) 'NEW  C(',jk,')  = ', ftempc * fthk 
    3234                            endif 
    3235 # endif 
    3236                         else 
    3237                         fq1      = fq0 * exp(-(fthk / xfastc))       !! how much organic C leaves this box        (mol) 
    3238                         freminc  = (fq0 - fq1) / fthk                !! C remineralisation in this box            (mol) 
    3239                         ffastc(ji,jj)  = fq1 
    3240                      endif 
    3241                      !! 
    3242                      !! === organic nitrogen === 
    3243                      fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol) 
    3244                      if (iball.eq.1) then 
    3245                         fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic N   (non-dim) 
    3246                         fq6      = (fq0 * fq5)                       !! how much organic N is unprotected         (mol) 
    3247                         fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected N leaves this box    (mol) 
    3248                         fq8      = (fq7 + (fq0 * fprotf))            !! how much total N leaves this box          (mol) 
    3249                         freminn  = (fq0 - fq8) / fthk                !! N remineralisation in this box            (mol) 
    3250                         ffastn(ji,jj)  = fq8 
    3251 # if defined key_debug_medusa 
    3252                         !! report in/out/remin fluxes of carbon for this level 
    3253                         if (idf.eq.1.AND.idfval.eq.1) then 
    3254                            IF (lwp) write (numout,*) '------------------------------' 
    3255                            IF (lwp) write (numout,*) 'totalN(',jk,')  = ', fq1 
    3256                            IF (lwp) write (numout,*) 'prtctN(',jk,')  = ', fq4 
    3257                            IF (lwp) write (numout,*) 'fprotf(',jk,')  = ', fprotf 
    3258                            IF (lwp) write (numout,*) '------------------------------' 
    3259                            if (freminn < 0.0) then 
    3260                               IF (lwp) write (numout,*) '** FREMIN ERROR **' 
    3261                            endif 
    3262                            IF (lwp) write (numout,*) 'IN   N(',jk,')  = ', fq0 
    3263                            IF (lwp) write (numout,*) 'LOST N(',jk,')  = ', freminn * fthk 
    3264                            IF (lwp) write (numout,*) 'OUT  N(',jk,')  = ', fq8 
    3265                            IF (lwp) write (numout,*) 'NEW  N(',jk,')  = ', ftempn * fthk 
    3266                         endif 
    3267 # endif 
    3268                      else 
    3269                         fq1      = fq0 * exp(-(fthk / xfastc))       !! how much organic N leaves this box        (mol) 
    3270                         freminn  = (fq0 - fq1) / fthk                !! N remineralisation in this box            (mol) 
    3271                         ffastn(ji,jj)  = fq1 
    3272                      endif 
    3273                      !! 
    3274                      !! === organic iron === 
    3275                      fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol) 
    3276                      if (iball.eq.1) then 
    3277                         fq5      = (1.0 - fprotf)                    !! unprotected fraction of total organic Fe  (non-dim) 
    3278                         fq6      = (fq0 * fq5)                       !! how much organic Fe is unprotected        (mol) 
    3279                         fq7      = (fq6 * exp(-(fthk / xfastc)))     !! how much unprotected Fe leaves this box   (mol) 
    3280                         fq8      = (fq7 + (fq0 * fprotf))            !! how much total Fe leaves this box         (mol)             
    3281                         freminfe = (fq0 - fq8) / fthk                !! Fe remineralisation in this box           (mol) 
    3282                         ffastfe(ji,jj) = fq8 
    3283                      else 
    3284                         fq1      = fq0 * exp(-(fthk / xfastc))       !! how much total Fe leaves this box         (mol)       
    3285                         freminfe = (fq0 - fq1) / fthk                !! Fe remineralisation in this box           (mol) 
    3286                         ffastfe(ji,jj) = fq1 
    3287                      endif 
    3288                      !! 
    3289                      !! === biogenic silicon === 
    3290                      fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol)  
    3291                      fq1      = fq0 * exp(-(fthk / xfastsi))         !! how much  opal leaves this box            (mol) 
    3292                      freminsi = (fq0 - fq1) / fthk                   !! Si remineralisation in this box           (mol) 
    3293                      ffastsi(ji,jj) = fq1 
    3294                      !! 
    3295                      !! === biogenic calcium carbonate === 
    3296                      fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol) 
    3297                      if (fdep.le.fccd_dep) then 
    3298                         !! whole grid cell above CCD 
    3299                         fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol) 
    3300                         freminca = 0.0                               !! above lysocline, no Ca dissolves          (mol) 
    3301                         fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#) 
    3302                      elseif (fdep.ge.fccd_dep) then 
    3303                         !! whole grid cell below CCD 
    3304                         fq1      = fq0 * exp(-(fthk / xfastca))      !! how much CaCO3 leaves this box            (mol) 
    3305                         freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol) 
    3306                      else 
    3307                         !! partial grid cell below CCD 
    3308                         fq2      = fdep1 - fccd_dep                  !! amount of grid cell below CCD             (m) 
    3309                         fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol) 
    3310                         freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol) 
    3311                      endif 
    3312                      ffastca(ji,jj) = fq1  
    3313                   else 
    3314                      !! this is BELOW THE LAST OCEAN BOX (do nothing) 
    3315                      freminc  = 0.0 
    3316                      freminn  = 0.0 
    3317                      freminfe = 0.0 
    3318                      freminsi = 0.0 
    3319                      freminca = 0.0               
    3320                   endif 
    3321  
    3322                elseif (jexport.eq.2.or.jexport.eq.3) then 
    3323                   if (jexport.eq.2) then 
    3324                      !!====================================================================== 
    3325                      !! MARTIN ET AL. (1987) SUBMODEL 
    3326                      !!====================================================================== 
    3327                      !!  
    3328                      !!---------------------------------------------------------------------- 
    3329                      !! This submodel uses the classic Martin et al. (1987) curve 
    3330                      !! to determine the attenuation of fast-sinking detritus down 
    3331                      !! the water column.  All three organic elements, C, N and Fe, 
    3332                      !! are handled identically, and their quantities in sinking 
    3333                      !! particles attenuate according to a power relationship 
    3334                      !! governed by parameter "b".  This is assigned a canonical  
    3335                      !! value of -0.858.  Biogenic opal and calcium carbonate are 
    3336                      !! attentuated using the same function as in the ballast 
    3337                      !! submodel 
    3338                      !!---------------------------------------------------------------------- 
    3339                      !! 
    3340                      fb_val = -0.858 
    3341                   elseif (jexport.eq.3) then 
    3342                      !!====================================================================== 
    3343                      !! HENSON ET AL. (2011) SUBMODEL 
    3344                      !!====================================================================== 
    3345                      !! 
    3346                      !!---------------------------------------------------------------------- 
    3347                      !! This submodel reconfigures the Martin et al. (1987) curve by 
    3348                      !! allowing the "b" value to vary geographically.  Its value is 
    3349                      !! set, following Henson et al. (2011), as a function of local 
    3350                      !! sea surface temperature: 
    3351                      !!   b = -1.06 + (0.024 * SST) 
    3352                      !! This means that remineralisation length scales are longer in 
    3353                      !! warm, tropical areas and shorter in cold, polar areas.  This 
    3354                      !! does seem back-to-front (i.e. one would expect GREATER 
    3355                      !! remineralisation in warmer waters), but is an outcome of  
    3356                      !! analysis of sediment trap data, and it may reflect details 
    3357                      !! of ecosystem structure that pertain to particle production 
    3358                      !! rather than simply Q10. 
    3359                      !!---------------------------------------------------------------------- 
    3360                      !! 
    3361                      fl_sst = tsn(ji,jj,1,jp_tem) 
    3362                      fb_val = -1.06 + (0.024 * fl_sst) 
    3363                   endif 
    3364                   !!    
    3365                   if (jk.eq.1) then 
    3366                      !! this is the SURFACE OCEAN BOX (no remineralisation) 
    3367                      !! 
    3368                      freminc  = 0.0 
    3369                      freminn  = 0.0 
    3370                      freminfe = 0.0 
    3371                      freminsi = 0.0 
    3372                      freminca = 0.0 
    3373                   elseif (jk.le.jmbathy) then 
    3374                      !! this is an OCEAN BOX (remineralise some material) 
    3375                      !! 
    3376                      !! === organic carbon === 
    3377                      fq0      = ffastc(ji,jj)                        !! how much organic C enters this box        (mol) 
    3378                      fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic C leaves this box        (mol) 
    3379                      freminc  = (fq0 - fq1) / fthk                   !! C remineralisation in this box            (mol) 
    3380                      ffastc(ji,jj)  = fq1 
    3381                      !! 
    3382                      !! === organic nitrogen === 
    3383                      fq0      = ffastn(ji,jj)                        !! how much organic N enters this box        (mol) 
    3384                      fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic N leaves this box        (mol) 
    3385                      freminn  = (fq0 - fq1) / fthk                   !! N remineralisation in this box            (mol) 
    3386                      ffastn(ji,jj)  = fq1 
    3387                      !! 
    3388                      !! === organic iron === 
    3389                      fq0      = ffastfe(ji,jj)                       !! how much organic Fe enters this box       (mol) 
    3390                      fq1      = fq0 * ((fdep1/fdep)**fb_val)         !! how much organic Fe leaves this box       (mol) 
    3391                      freminfe = (fq0 - fq1) / fthk                   !! Fe remineralisation in this box           (mol) 
    3392                      ffastfe(ji,jj) = fq1 
    3393                      !! 
    3394                      !! === biogenic silicon === 
    3395                      fq0      = ffastsi(ji,jj)                       !! how much  opal centers this box           (mol)  
    3396                      fq1      = fq0 * exp(-(fthk / xfastsi))         !! how much  opal leaves this box            (mol) 
    3397                      freminsi = (fq0 - fq1) / fthk                   !! Si remineralisation in this box           (mol) 
    3398                      ffastsi(ji,jj) = fq1 
    3399                      !! 
    3400                      !! === biogenic calcium carbonate === 
    3401                      fq0      = ffastca(ji,jj)                       !! how much CaCO3 enters this box            (mol) 
    3402                      if (fdep.le.ocal_ccd(ji,jj)) then 
    3403                         !! whole grid cell above CCD 
    3404                         fq1      = fq0                               !! above lysocline, no Ca dissolves          (mol) 
    3405                         freminca = 0.0                               !! above lysocline, no Ca dissolves          (mol) 
    3406                         fccd(ji,jj) = real(jk)                       !! which is the last level above the CCD?    (#) 
    3407                      elseif (fdep.ge.ocal_ccd(ji,jj)) then 
    3408                         !! whole grid cell below CCD 
    3409                         fq1      = fq0 * exp(-(fthk / xfastca))      !! how much CaCO3 leaves this box            (mol) 
    3410                         freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol) 
    3411                      else 
    3412                         !! partial grid cell below CCD 
    3413                         fq2      = fdep1 - ocal_ccd(ji,jj)           !! amount of grid cell below CCD             (m) 
    3414                         fq1      = fq0 * exp(-(fq2 / xfastca))       !! how much CaCO3 leaves this box            (mol) 
    3415                         freminca = (fq0 - fq1) / fthk                !! Ca remineralisation in this box           (mol) 
    3416                      endif 
    3417                      ffastca(ji,jj) = fq1  
    3418                   else 
    3419                      !! this is BELOW THE LAST OCEAN BOX (do nothing) 
    3420                      freminc  = 0.0 
    3421                      freminn  = 0.0 
    3422                      freminfe = 0.0 
    3423                      freminsi = 0.0 
    3424                      freminca = 0.0               
    3425                   endif 
    3426  
    3427                endif 
    3428  
    3429                !!---------------------------------------------------------------------- 
    3430                !! Fast-sinking detritus fluxes, pt. 2: UPDATE FAST FLUXES 
    3431                !! here locally calculated additions to the fast-sinking flux are added 
    3432                !! to the total fast-sinking flux; this is done here such that material 
    3433                !! produced in a particular layer is only remineralised below this 
    3434                !! layer 
    3435                !!---------------------------------------------------------------------- 
    3436                !! 
    3437                !! add sinking material generated in this layer to running totals 
    3438                !! 
    3439                !! === organic carbon ===                          (diatom and mesozooplankton mortality) 
    3440                ffastc(ji,jj)  = ffastc(ji,jj)  + (ftempc  * fthk) 
    3441                !! 
    3442                !! === organic nitrogen ===                        (diatom and mesozooplankton mortality) 
    3443                ffastn(ji,jj)  = ffastn(ji,jj)  + (ftempn  * fthk) 
    3444                !! 
    3445                !! === organic iron ===                            (diatom and mesozooplankton mortality) 
    3446                ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe * fthk) 
    3447                !! 
    3448                !! === biogenic silicon ===                        (diatom mortality and grazed diatoms) 
    3449                ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi * fthk) 
    3450                !! 
    3451                !! === biogenic calcium carbonate ===              (latitudinally-based fraction of total primary production) 
    3452                ffastca(ji,jj) = ffastca(ji,jj) + (ftempca * fthk) 
    3453  
    3454                !!---------------------------------------------------------------------- 
    3455                !! Fast-sinking detritus fluxes, pt. 3: SEAFLOOR 
    3456                !! remineralise all remaining fast-sinking detritus to dissolved 
    3457                !! nutrients; the sedimentation fluxes calculated here allow the 
    3458                !! separation of what's remineralised sinking through the final 
    3459                !! ocean box from that which is added to the final box by the 
    3460                !! remineralisation of material that reaches the seafloor (i.e. 
    3461                !! the model assumes that *all* material that hits the seafloor 
    3462                !! is remineralised and that none is permanently buried; hey, 
    3463                !! this is a giant GCM model that can't be run for long enough 
    3464                !! to deal with burial fluxes!) 
    3465                !! 
    3466                !! in a change to this process, in part so that MEDUSA behaves 
    3467                !! a little more like ERSEM et al., fast-sinking detritus (N, Fe 
    3468                !! and C) is converted to slow sinking detritus at the seafloor 
    3469                !! instead of being remineralised; the rationale is that in 
    3470                !! shallower shelf regions (... that are not fully mixed!) this 
    3471                !! allows the detrital material to return slowly to dissolved  
    3472                !! nutrient rather than instantaneously as now; the alternative 
    3473                !! would be to explicitly handle seafloor organic material - a 
    3474                !! headache I don't wish to experience at this point; note that 
    3475                !! fast-sinking Si and Ca detritus is just remineralised as  
    3476                !! per usual 
    3477                !!  
    3478                !! AXY (13/01/12) 
    3479                !! in a further change to this process, again so that MEDUSA is 
    3480                !! a little more like ERSEM et al., material that reaches the 
    3481                !! seafloor can now be added to sediment pools and stored for 
    3482                !! slow release; there are new 2D arrays for organic nitrogen, 
    3483                !! iron and carbon and inorganic silicon and carbon that allow 
    3484                !! fast and slow detritus that reaches the seafloor to be held 
    3485                !! and released back to the water column more slowly; these arrays 
    3486                !! are transferred via the tracer restart files between repeat 
    3487                !! submissions of the model 
    3488                !!---------------------------------------------------------------------- 
    3489                !!  
    3490                ffast2slowc  = 0.0 
    3491                ffast2slown  = 0.0 
    3492                ffast2slowfe = 0.0 
    3493                !! 
    3494                if (jk.eq.jmbathy) then 
    3495                   !! this is the BOTTOM OCEAN BOX (remineralise everything) 
    3496                   !! 
    3497                   !! AXY (17/01/12): tweaked to include benthos pools 
    3498                   !!  
    3499                   !! === organic carbon === 
    3500                   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
    3501                      freminc  = freminc + (ffastc(ji,jj) / fthk)    !! C remineralisation in this box            (mol/m3) 
    3502                   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
    3503                      ffast2slowc = ffastc(ji,jj) / fthk             !! fast C -> slow C                          (mol/m3) 
    3504                      fslowc      = fslowc + ffast2slowc 
    3505                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
    3506                      f_fbenin_c(ji,jj)  = ffastc(ji,jj)             !! fast C -> benthic C                       (mol/m2) 
    3507                   endif 
    3508                   fsedc(ji,jj)   = ffastc(ji,jj)                          !! record seafloor C                         (mol/m2) 
    3509                   ffastc(ji,jj)  = 0.0 
    3510                   !! 
    3511                   !! === organic nitrogen === 
    3512                   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
    3513                      freminn  = freminn + (ffastn(ji,jj) / fthk)    !! N remineralisation in this box            (mol/m3) 
    3514                   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
    3515                      ffast2slown = ffastn(ji,jj) / fthk             !! fast N -> slow N                          (mol/m3) 
    3516                      fslown      = fslown + ffast2slown 
    3517                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
    3518                      f_fbenin_n(ji,jj)  = ffastn(ji,jj)             !! fast N -> benthic N                       (mol/m2) 
    3519                   endif 
    3520                   fsedn(ji,jj)   = ffastn(ji,jj)                          !! record seafloor N                         (mol/m2) 
    3521                   ffastn(ji,jj)  = 0.0 
    3522                   !! 
    3523                   !! === organic iron === 
    3524                   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
    3525                      freminfe = freminfe + (ffastfe(ji,jj) / fthk)  !! Fe remineralisation in this box           (mol/m3) 
    3526                   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
    3527                      ffast2slowfe = ffastn(ji,jj) / fthk            !! fast Fe -> slow Fe                        (mol/m3) 
    3528                   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
    3529                      f_fbenin_fe(ji,jj) = ffastfe(ji,jj)            !! fast Fe -> benthic Fe                     (mol/m2) 
    3530                   endif 
    3531                   fsedfe(ji,jj)  = ffastfe(ji,jj)                         !! record seafloor Fe                        (mol/m2) 
    3532                   ffastfe(ji,jj) = 0.0 
    3533                   !! 
    3534                   !! === biogenic silicon === 
    3535                   if (jinorgben.eq.0) then 
    3536                      freminsi = freminsi + (ffastsi(ji,jj) / fthk)  !! Si remineralisation in this box           (mol/m3) 
    3537                   elseif (jinorgben.eq.1) then 
    3538                      f_fbenin_si(ji,jj) = ffastsi(ji,jj)            !! fast Si -> benthic Si                     (mol/m2) 
    3539                   endif 
    3540                   fsedsi(ji,jj)   = ffastsi(ji,jj)                         !! record seafloor Si                        (mol/m2) 
    3541                   ffastsi(ji,jj) = 0.0 
    3542                   !! 
    3543                   !! === biogenic calcium carbonate === 
    3544                   if (jinorgben.eq.0) then 
    3545                      freminca = freminca + (ffastca(ji,jj) / fthk)  !! Ca remineralisation in this box           (mol/m3)  
    3546                   elseif (jinorgben.eq.1) then 
    3547                      f_fbenin_ca(ji,jj) = ffastca(ji,jj)            !! fast Ca -> benthic Ca                     (mol/m2) 
    3548                   endif 
    3549                   fsedca(ji,jj)   = ffastca(ji,jj)                         !! record seafloor Ca                        (mol/m2) 
    3550                   ffastca(ji,jj) = 0.0 
    3551                endif 
    3552  
    3553 # if defined key_debug_medusa 
    3554                if (idf.eq.1) then 
    3555                   !!---------------------------------------------------------------------- 
    3556                   !! Integrate total fast detritus remineralisation 
    3557                   !!---------------------------------------------------------------------- 
    3558                   !! 
    3559                   fofd_n(ji,jj)  = fofd_n(ji,jj)  + (freminn  * fthk) 
    3560                   fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi * fthk) 
    3561                   fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe * fthk) 
    3562 #  if defined key_roam 
    3563                   fofd_c(ji,jj)  = fofd_c(ji,jj)  + (freminc  * fthk) 
    3564 #  endif 
    3565                endif 
    3566 # endif 
    3567  
    3568                !!---------------------------------------------------------------------- 
    3569                !! Sort out remineralisation tally of fast-sinking detritus 
    3570                !!---------------------------------------------------------------------- 
    3571                !! 
    3572                !! update fast-sinking regeneration arrays 
    3573                fregenfast(ji,jj)   = fregenfast(ji,jj)   + (freminn  * fthk) 
    3574                fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi * fthk) 
    3575 # if defined key_roam 
    3576                fregenfastc(ji,jj)  = fregenfastc(ji,jj)  + (freminc  * fthk) 
    3577 # endif 
    3578  
    3579                !!---------------------------------------------------------------------- 
    3580                !! Benthic remineralisation fluxes 
    3581                !!---------------------------------------------------------------------- 
    3582                !! 
    3583                if (jk.eq.jmbathy) then 
    3584                   !! 
    3585                   !! organic components 
    3586                   if (jorgben.eq.1) then 
    3587                      f_benout_n(ji,jj)  = xsedn  * zn_sed_n(ji,jj) 
    3588                      f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj) 
    3589                      f_benout_c(ji,jj)  = xsedc  * zn_sed_c(ji,jj) 
    3590                   endif 
    3591                   !! 
    3592                   !! inorganic components 
    3593                   if (jinorgben.eq.1) then 
    3594                      f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj) 
    3595                      f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) 
    3596                      !! 
    3597                      !! account for CaCO3 that dissolves when it shouldn't 
    3598                      if ( fdep .le. fccd_dep ) then 
    3599                         f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) 
    3600                      endif 
    3601                   endif 
    3602                endif 
    3603                CALL flush(numout) 
    3604  
    3605                !!====================================================================== 
    3606                !! LOCAL GRID CELL TRENDS 
    3607                !!====================================================================== 
    3608                !! 
    3609                !!---------------------------------------------------------------------- 
    3610                !! Determination of trends 
    3611                !!---------------------------------------------------------------------- 
    3612                !! 
    3613                !!---------------------------------------------------------------------- 
    3614                !! chlorophyll 
    3615                btra(jpchn) = b0 * ( & 
    3616                  + ((frn * fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2) * (fthetan / xxi) ) 
    3617                btra(jpchd) = b0 * ( & 
    3618                  + ((frd * fprd * zphd) - fgmepd - fdpd - fdpd2) * (fthetad / xxi) ) 
    3619                !! 
    3620                !!---------------------------------------------------------------------- 
    3621                !! phytoplankton 
    3622                btra(jpphn) = b0 * ( & 
    3623                  + (fprn * zphn) - fgmipn - fgmepn - fdpn - fdpn2 ) 
    3624                btra(jpphd) = b0 * ( & 
    3625                  + (fprd * zphd) - fgmepd - fdpd - fdpd2 ) 
    3626                btra(jppds) = b0 * ( & 
    3627                  + (fprds * zpds) - fgmepds - fdpds - fsdiss - fdpds2 ) 
    3628                !! 
    3629                !!---------------------------------------------------------------------- 
    3630                !! zooplankton 
    3631                btra(jpzmi) = b0 * ( & 
    3632                  + fmigrow - fgmezmi - fdzmi - fdzmi2 ) 
    3633                btra(jpzme) = b0 * ( & 
    3634                  + fmegrow - fdzme - fdzme2 ) 
    3635                !! 
    3636                !!---------------------------------------------------------------------- 
    3637                !! detritus 
    3638                btra(jpdet) = b0 * ( & 
    3639                  + fdpn + ((1.0 - xfdfrac1) * fdpd)              &  ! mort. losses 
    3640                  + fdzmi + ((1.0 - xfdfrac2) * fdzme)            &  ! mort. losses 
    3641                  + ((1.0 - xbetan) * (finmi + finme))            &  ! assim. inefficiency 
    3642                  - fgmid - fgmed - fdd                           &  ! grazing and remin. 
    3643                  + ffast2slown )