SUBROUTINE ice_bio_diff(kideb,kiut,nlay_i) !------------------------------------------------------------------------------! ! ! --- ice_bio_diff --- ! ! Transport and diffusion of tracers ! (c) Martin Vancoppenolle, May 2007 ! 1.1 Rayleigh number based diffusivity, Oct 2008 ! 1.2 Includes open upper boundary condition for gas exchange Jul 2010 - Mar 2011 ! 1.3 Simplification and implementation of flooding velocity, Feb 2012 ! ! Rewriting of surface gas boundary condition (M. Vancoppenolle & S. Moreau, 2015) ! - routine for gas transfer velocity (ice_gas_trvel) ! - add solubility and partial pressure for each gas ! ! ! Use the real surface pCO2 and not the 1st layer value ! ! remove ji ! !------------------------------------------------------------------------------! USE lib_fortran INCLUDE 'type.com' INCLUDE 'para.com' INCLUDE 'const.com' INCLUDE 'ice.com' INCLUDE 'thermo.com' INCLUDE 'bio.com' INTEGER :: & ji , ! : horizontal space index & jn ! : horizontal space index jn REAL(8), DIMENSION( maxnlay ) :: !: dummy factors for tracer equation & za , !: winter & zb , & zc , & ze , !: summer & zind , !: independent term in the tridiag system & zindtbis , !: & zdiagbis REAL(8), DIMENSION(maxnlay,3) ::!: dummy factors for tracer equation & ztrid !: tridiagonal matrix REAL(8), DIMENSION(nlay_bio) :: & z_sbr_i !: brine salinity REAL(8), DIMENSION(ntra_bio_max) :: & zpp_gas !: partial pressure REAL(8) :: & zdummy1 , !: dummy factors & zdummy2 , !: & zdummy3 , !: & zswitchs , !: switch for summer drainage & f_sn_rat , & wspd_trs , & zsat_arg , & zsat_oxy , & zsat_CO2 , & zsat_nit , & mol_diff INTEGER :: & indtr , !: index of tridiagonal system & iter !: time step index LOGICAL :: & ln_write_bio , & ln_con_bio , & ln_flood ln_write_bio = .TRUE. ln_con_bio = .TRUE. ln_flood = .TRUE. zsol = 0. zb0 = 0. zpatm_gas = 0. zpp_gas(:) = 0.0 !======================================================================= WRITE(numout,*) WRITE(numout,*) ' ** ice_bio_diff : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) DO ji = kideb, kiut ! !----------------------------------------------------------------------- ! 1) Initialization !----------------------------------------------------------------------- ! IF ( ln_write_bio ) THEN WRITE(numout,*) ' Initialization ' WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' WRITE(numout,*) ENDIF !--------------- ! Interpolation !--------------- CALL ice_bio_interp_phy2bio(kideb,kiut,nlay_i,.FALSE.) ! interpolation of physical variables ! on the biological grid ! mass of salt, heat content, brine volume, Rb, PAR CALL ice_bio_interp_diffus(kideb,kiut,nlay_i,.TRUE.) IF ( ln_write_bio ) THEN WRITE(numout,*) ' diff_br_bio : ', ( diff_br_bio(layer), & layer = 1, nlay_bio ) ENDIF !-------------------------------- ! Brine concentration of tracers !-------------------------------- DO jn = 1, ntra_bio IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN DO jk = 1, nlay_bio c_i_bio(jn,jk) = cbu_i_bio(jn,jk) / e_i_bio(jk) END DO ENDIF IF ( flag_adsorb(jn) .AND. flag_active(jn) ) THEN DO jk = 1, nlay_bio c_i_bio(jn,jk) = 0. END DO ENDIF WRITE(numout,*) ' 01 *** jn = ', jn WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:) END DO !--------------------------------------------------------------- ! Equilibrate carbonate system to get aqueous CO2 concentration !--------------------------------------------------------------- IF ( ln_carbon ) CALL ice_carb_chem DO jn = 1, ntra_bio WRITE(numout,*) ' 02 *** jn = ', jn WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:) END DO !-------------------- ! Conservation check !-------------------- CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_init,cbu_i_bio, & deltaz_i_bio, .FALSE.) IF ( ln_write_bio ) THEN DO jn = 1, ntra_bio WRITE(numout,*) ' mt_i_bio_init : ', mt_i_bio_init(jn) END DO ENDIF ! layer by layer DO jn = 1, ntra_bio DO jk = 1, nlay_bio m_i_bio_init(jn,jk) = cbu_i_bio(jn,jk)*deltaz_i_bio(jk) END DO END DO ! !----------------------------------------------------------------------- ! 3) Surface boundary condition for gases !----------------------------------------------------------------------- ! zaeff = e_i_bio(1) * brines_ar ! proportionality factor for surface flux zpp_gas(:) = mixr_gas(:) * psbqb(ji) / ! atmospheric partial pressure & 101325. !---------------------------- CALL ice_gas_solu ! Solubilities (sol_gas(jn)) !---------------------------- !------------------------- gas_trvel(:) = 0. ! Gas transfer velocities !------------------------- DO jn = 1, ntra_bio IF ( flag_active(jn) .AND. ( biotr_i_typ(jn) .EQ. 'gas' ) .AND. ! select gases only & ( flag_diff(jn) .OR. flag_adsorb(jn) ) .AND. ln_gasflux ) THEN SELECT CASE ( i_gasflux ) CASE (1) ! No flux CASE (2) ! Gas diffusion between ice and atm IF ( e_i_bio(1) .GT. e_thr_gasflux ) & gas_trvel(jn) = dmol_gas(jn) / h_bl_gas END SELECT IF ( ln_write_bio ) THEN WRITE(numout,*) ' *** Gas transfer velocity --- ' WRITE(numout,*) ' --- Tracer --- : ', biotr_i_nam(jn) WRITE(numout,*) ' gas_trvel : ', gas_trvel(jn) ENDIF ! ln_write_bio ENDIF ! flags ! !----------------------------------------------------------------------- ! 3) Switches, flushing and flooding velocities !----------------------------------------------------------------------- ! !---------- ! Switches !---------- ! summer switch zbvmin = 1.0 DO layer = 1, nlay_bio zbvmin = MIN( e_i_bio(layer) , zbvmin ) ! minimum brine volume END DO zswitchs = MAX( 0.0, SIGN ( 1.0d0, t_su_b(ji) - tpw ) ) ! 0 si hiver 1 si ete IF ( zbvmin .LT. e_thr_flu ) zswitchs = 0.0 IF ( ln_write_bio ) THEN WRITE(numout,*) ' zswitchs : ', zswitchs WRITE(numout,*) ENDIF ! flood switch IF ( w_flood .GT. 0 ) THEN z_flood = 0.0 ELSE z_flood = 1.0 ENDIF ! recompute flushing velocity w_flush = qsummer / ddtb / e_i_bio(1) ! !----------------------------------------------------------------------- ! 4) Compute dummy factors for tracer diffusion equation !----------------------------------------------------------------------- ! IF ( ( flag_diff(jn) .OR. flag_adsorb(jn) ) & .AND. flag_active(jn) ) THEN IF ( ln_write_bio ) THEN WRITE(numout,*) ' --------------------------------- ' WRITE(numout,*) ' Diffusion for ', biotr_i_nam(jn) WRITE(numout,*) ' --------------------------------- ' WRITE(numout,*) WRITE(numout,*) ' Dummy factors ' WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' WRITE(numout,*) ENDIF !-------------------- ! za factors !-------------------- DO layer = 1, nlay_bio za(layer) = ddtb / ( deltaz_i_bio(layer) * e_i_bio(layer) ) END DO !-------------------- ! zb, zc, ze factors !-------------------- DO layer = 1, nlay_bio - 1 ! interpolate brine volume at the interface between layers zdummy1 = ( e_i_bio(layer + 1) - e_i_bio(layer) ) / & ( z_i_bio(layer + 1) - z_i_bio(layer) ) zdummy2 = deltaz_i_bio(layer) / 2.0 zdummy3 = e_i_bio(layer) + zdummy1 * zdummy2 zb(layer) = zdummy3 * diff_br_bio(layer) / & ( z_i_bio(layer+1) - z_i_bio(layer) ) zc(layer) = w_flood * zdummy3 * z_flood ze(layer) = ( w_flood * ( 1. - z_flood ) + w_flush ) * zdummy3 END DO ! Basal values zb(nlay_bio) = 2. * e_i_bio(nlay_bio) / & deltaz_i_bio(nlay_bio) * diff_br_bio(nlay_bio) zc(nlay_bio) = w_flood * e_i_bio(nlay_bio) * z_flood ze(nlay_bio) = ( w_flood * ( 1. - z_flood ) + w_flush ) * & e_i_bio(nlay_bio) ! Open upper boundary condition for gas zb0 = 0. IF ( biotr_i_typ(jn) .EQ. 'gas' ) zb0 = gas_trvel(jn) * zaeff ! Block fluxes above the biologically active layer (SL & BAL cases) IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) THEN zb(1:nlay_bio-1) = 0. zc(1:nlay_bio-1) = 0. ze(1:nlay_bio-1) = 0. ENDIF IF ( ln_write_bio ) THEN WRITE(numout,*) ' Winter factors ' WRITE(numout,*) ' za : ', ( za (layer), & layer = 1, nlay_bio) WRITE(numout,*) ' zb : ', ( zb (layer), & layer = 1, nlay_bio) WRITE(numout,*) ' zc : ', ( zc (layer), & layer = 1, nlay_bio) WRITE(numout,*) ' ze : ', ( ze (layer), & layer = 1, nlay_bio) WRITE(numout,*) ENDIF ! !----------------------------------------------------------------------- ! 5) Tridiagonal system terms for tracer diffusion equation, winter !----------------------------------------------------------------------- ! !---------------- ! first equation !---------------- ztrid(1,1) = 0.0 ztrid(1,2) = 1.0 + za(1) * ( zb(1) + ze(1) + zb0 ) ztrid(1,3) = za(1) * ( -zb(1) + zc(1) ) zind(1) = c_i_bio(jn,1) + za(1)*zb0*sol_gas(jn,1)*zpp_gas(jn) ! IF ( jn .NE. jn_dic ) THEN ! zind(1) = c_i_bio(jn,1) + za(1)*zb0*sol_gas(jn,1)*zpp_gas(jn) ! ELSE ! CALL ice_carb_chem ! zind(1) = c_i_bio(jn,1) + za(1)*zb0*( sol_gas(jn,1)*zpp_gas(jn) - co2aq(1) ) ! ENDIF ! for dic would read c_i_bio(jn,1) + za(1)*zb0*(sol_gas(jn,1)*zpp_gas(jn) - zco2 ) ! where zco2 = co2aq(1) !----------------- ! inner equations !----------------- DO layer = 2, nlay_bio - 1 ztrid(layer,1) = - za(layer) * ( zb(layer-1) + ze(layer-1) ) ztrid(layer,2) = 1.0 + za(layer) * ( zb(layer) + & ze(layer) + zb(layer-1) - zc(layer-1) ) ztrid(layer,3) = za(layer) * ( -zb(layer) + zc(layer) ) zind(layer) = c_i_bio(jn,layer) END DO !---------------- ! last equation !---------------- ztrid(nlay_bio,1) = -za(nlay_bio) * ( zb(nlay_bio-1) + & ze(nlay_bio-1) ) ztrid(nlay_bio,2) = 1.0 + za(nlay_bio) * ( zb(nlay_bio) + & ze(nlay_bio) + zb(nlay_bio-1) - & zc(nlay_bio-1) ) ztrid(nlay_bio,3) = 0. zind(nlay_bio) = c_i_bio(jn,nlay_bio) + za(nlay_bio) * & ( zb(nlay_bio) - zc(nlay_bio) )*c_w_bio(jn) IF ( ln_write_bio ) THEN WRITE(numout,*) WRITE(numout,*) ' Tridiag terms : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~ ' WRITE(numout,*) DO layer = 1, nlay_bio WRITE(numout,*) ' layer : ', layer WRITE(numout,*) ' ztrid : ', ztrid(layer,1), & ztrid(layer,2), & ztrid(layer,3) WRITE(numout,*) ' zind : ',zind(layer) END DO ENDIF ! !----------------------------------------------------------------------- ! 6) Solving the tridiagonal system !----------------------------------------------------------------------- ! ! The tridiagonal system is solved with Gauss elimination ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, ! McGraw-Hill 1984. zindtbis(1) = zind(1) zdiagbis(1) = ztrid(1,2) DO layer = 2, nlay_bio zdiagbis(layer) = ztrid(layer,2) - ztrid(layer,1) * & ztrid(layer-1,3) / zdiagbis(layer-1) zindtbis(layer) = zind(layer) - ztrid(layer,1) * & zindtbis(layer-1) / zdiagbis(layer-1) END DO !----------------------------- ! Tracer brine concentrations !----------------------------- c_i_bio(jn,nlay_bio) = zindtbis(nlay_bio) / zdiagbis(nlay_bio) DO layer = nlay_bio - 1 , 1 , -1 c_i_bio(jn,layer) = (zindtbis(layer) - ztrid(layer,3)* & c_i_bio(jn,layer+1)) / zdiagbis(layer) END DO IF ( ln_write_bio ) THEN WRITE(numout,*) WRITE(numout,*) ' Resolution ' WRITE(numout,*) ' ~~~~~~~~~~~ ' WRITE(numout,*) DO layer = 1, nlay_bio WRITE(numout,*) ' layer : ', layer WRITE(numout,*) ' zdiagbis : ', zdiagbis(layer) WRITE(numout,*) ' zindtbis : ', zindtbis(layer) WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,layer) END DO ENDIF ENDIF ! flag ! Update gas flux fgas(jn) = 0. IF ( ( biotr_i_typ(jn) .EQ. 'gas' ) .AND. & flag_active(jn) ) THEN fgas(jn) = - zb0 * ( c_i_bio(jn,1) - & sol_gas(jn,1)*zpp_gas(jn) ) ! mmol/m2/s, positive to the ice ENDIF WRITE(numout,*) WRITE(numout,*) ' tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' fgas : ', fgas(jn) WRITE(numout,*) END DO ! jn ! !----------------------------------------------------------------------- ! 7) Recover bulk tracer concentrations !----------------------------------------------------------------------- ! !-------------------------------- ! Tracer bulk ice concentrations !-------------------------------- IF ( ln_write_bio ) THEN WRITE(numout,*) WRITE(numout,*) ' fgas Argon : ', fgas(jn_arg) WRITE(numout,*) ' fgas Oxygen : ', fgas(jn_oxy) WRITE(numout,*) ' fgas CO2 : ', fgas(jn_co2) WRITE(numout,*) ENDIF DO jn = 1, ntra_bio IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN DO layer = 1, nlay_bio cbu_i_bio(jn,layer) = c_i_bio(jn,layer) * e_i_bio(layer) END DO WRITE(numout,*) ' 03 *** jn = ', jn WRITE(numout,*) ' cbu_i_bio : ', cbu_i_bio(jn,:) WRITE(numout,*) ' c_i_bio : ', c_i_bio(jn,:) IF ( jn .EQ. jn_dic ) THEN WRITE(numout,*) WRITE(numout,*) ' Change in DIC :' WRITE(numout,*) ' c_i_bio(11,1) before :',c_i_bio(jn,1) WRITE(numout,*) ' cbu_i_bio(11,1) befe :',cbu_i_bio(jn,1) cbu_i_bio(jn,1) = cbu_i_bio(jn,1) + fgas(jn_co2) * ! change in bulk DIC due to fgas & ddtb / deltaz_i_bio(1) c_i_bio(jn,1) = cbu_i_bio(jn,1) / e_i_bio(1) ! update surface brine DIC after CO2 flux WRITE(numout,*) ' c_i_bio(11,1) after :',c_i_bio(jn,1) WRITE(numout,*) ' cbu_i_bio(11,1) after:',cbu_i_bio(jn,1) WRITE(numout,*) ' fgas(20) :', fgas(jn_co2) WRITE(numout,*) ENDIF ENDIF ! flag_diff IF ( flag_adsorb(jn) .AND. flag_active(jn) ) THEN DO layer = 1, nlay_bio cbu_i_bio(jn,layer) = cbu_i_bio(jn,layer) & + c_i_bio(jn,layer) * e_i_bio(layer) END DO ENDIF ! flag_adsorb END DO ! jn IF ( ln_write_bio ) THEN DO jn = 1, ntra_bio IF ( ( flag_diff(jn) .OR. flag_adsorb(jn) ) & .AND. flag_active(jn) ) THEN WRITE(numout,*) WRITE(numout,*) ' Tracer concentrations ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), & layer = 1, nlay_bio ) WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn,layer), & layer = 1, nlay_bio ) WRITE(numout,*) WRITE(numout,*) ' The ', biotr_i_nam(jn), 'fluxdwn is :', & fgas(jn), ' mol / m2 / s ' ENDIF ! flag_diff END DO ! jn ENDIF ! ln_write_bio ! !----------------------------------------------------------------------- ! 8) Conservation check !----------------------------------------------------------------------- ! IF ( ln_con_bio ) THEN IF ( ln_write_bio ) THEN WRITE(numout,*) WRITE(numout,*) ' Conservation check : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) ENDIF ! ln_write_bio CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_final,cbu_i_bio, & deltaz_i_bio, .FALSE.) zerror = 1.0d-15 DO jn = 1, ntra_bio IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN IF ( ln_write_bio ) THEN WRITE(numout,*) ' mt_i_bio_final : ', mt_i_bio_final(jn) ENDIF zfb = - e_i_bio(nlay_bio) * ( diff_br_bio(nlay_bio) * 2.0 / & deltaz_i_bio(nlay_bio) * ( c_i_bio(jn,nlay_bio) - & c_w_bio(jn) ) + w_flood * ( z_flood * c_w_bio(jn) + & ( 1. - z_flood ) * c_i_bio(jn,nlay_bio) ) + w_flush * & c_i_bio(jn,nlay_bio) ) f_bo_tra(jn) = zfb fcb(jn) = - zfb ! ice-ocean tracer flux f_su_tra(jn) = fgas(jn) & + zswitchs * ( qsummer * c_s_bio(jn) ) / ddtb ! fcb_max = idealized flux if c_i_bio(nlay_bio) .EQ. 0.0 all the time fcb_max(jn)= e_i_bio(nlay_bio) * ( diff_br_bio(nlay_bio) * 2.0 / & deltaz_i_bio(nlay_bio) * ( 0.0 - & c_w_bio(jn) ) + w_flood * ( z_flood * c_w_bio(jn) + & ( 1. - z_flood ) * 0. ) + w_flush * 0. ) IF ( ln_write_bio ) THEN WRITE(numout,*) ' f_bo_tra : ', f_bo_tra(jn) WRITE(numout,*) ' f_su_tra : ', f_su_tra(jn) WRITE(numout,*) ' zswitchs : ', zswitchs WRITE(numout,*) ' qsummer : ', qsummer WRITE(numout,*) ' c_i_bio(nlay_bio) : ',c_i_bio(jn,nlay_bio) WRITE(numout,*) ' mt_i_bio_init : ', mt_i_bio_init(jn) WRITE(numout,*) ' mt_i_bio_final : ', mt_i_bio_final(jn) WRITE(numout,*) ENDIF ! ln_write_bio ENDIF ! flag_diff END DO ! jn CALL ice_bio_conserv(kideb,kiut,ntra_bio,'ice_bio_diff ', & biotr_i_nam, zerror, & mt_i_bio_init,mt_i_bio_final, & f_bo_tra, f_su_tra, ddtb) ENDIF ! ln_con_bio ! !----------------------------------------------------------------------- ! X) Control Print !----------------------------------------------------------------------- ! WRITE(numout,*) WRITE(numout,*) ' ** After CO2 fluxes with atmosphere : ' WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' WRITE(numout,*) DO jn = 1, ntra_bio IF ( biotr_i_nam(jn) .EQ. 'CO2' ) THEN WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3) & layer = 1, nlay_bio ) WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice & layer = 1, nlay_bio ) ENDIF IF ( biotr_i_nam(jn) .EQ. 'DIC' ) THEN WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3) & layer = 1, nlay_bio ) WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice & layer = 1, nlay_bio ) ENDIF IF ( biotr_i_nam(jn) .EQ. 'Alk' ) THEN WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3) & layer = 1, nlay_bio ) WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice & layer = 1, nlay_bio ) ENDIF IF ( biotr_i_nam(jn) .EQ. 'Ika' ) THEN WRITE(numout,*) ' Tracer : ', biotr_i_nam(jn) WRITE(numout,*) ' c_i_bio : ', ( c_i_bio(jn,layer), ! concentration of tracers in brines (mmol m-3) & layer = 1, nlay_bio ) WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer), ! concentration of tracers in bulk ice & layer = 1, nlay_bio ) ENDIF END DO ! !----------------------------------------------------------------------- ! 8) Flux divergence diagnostics !----------------------------------------------------------------------- ! ! not sure if those are good anymore WRITE(numout,*) WRITE(numout,*) ' Diagnostics ' WRITE(numout,*) WRITE(numout,*) ! Right now, does not include summer DO jn = 1, ntra_bio IF ( flag_diff(jn) .AND. flag_active(jn) ) THEN fdiff(jn,0) = 0.0 ! top flux DO layer = 1, nlay_bio - 1 ! inner fluxes ! interpolate brine volume at the interface between layers zdummy1 = ( e_i_bio(layer + 1 ) - e_i_bio(layer) ) / & ( z_i_bio(layer + 1) - z_i_bio(layer) ) zdummy2 = deltaz_i_bio(layer) / 2.0 zdummy3 = e_i_bio(layer) + zdummy1 * zdummy2 fdiff(jn,layer) = - zdummy3 * & diff_br_bio(layer) * & ( c_i_bio(jn,layer+1) - c_i_bio(jn,layer) ) / & ( z_i_bio(layer + 1) - z_i_bio(layer) ) END DO ! lower flux fdiff(jn,nlay_bio) = - 2.0 * e_i_bio(nlay_bio) * & diff_br_bio(layer) * & ( c_w_bio(jn) - c_i_bio(jn,nlay_bio) ) / & deltaz_i_bio(nlay_bio) WRITE(numout,*) ' c_w_bio : ', c_w_bio(jn) WRITE(numout,*) ' c_i_bio(N) : ', c_i_bio(jn,nlay_bio) WRITE(numout,*) ' deltaz_i_bio : ', deltaz_i_bio(nlay_bio) ! divergence of the fluxes DO layer = 1, nlay_bio diag_divf_bio(jn,layer) = & - ( fdiff(jn,layer) - fdiff(jn,layer-1) ) / & deltaz_i_bio(layer) END DO DO layer = 1, nlay_bio - 1 WRITE(numout,*) ' fdiff(layer) : ', fdiff(jn,layer) END DO DO layer = 1, nlay_bio - 1 WRITE(numout,*) ' div_F(layer) : ', diag_divf_bio(jn,layer) END DO ENDIF ! flag_diff END DO ! jn ! !----------------------------------------------------------------------- ! 9) End of the routine !----------------------------------------------------------------------- ! END DO ! ji IF ( ln_write_bio ) THEN WRITE(numout,*) WRITE(numout,*) ' *** After diffusion of tracers *** ' WRITE(numout,*) ' model output ' DO jn = 1, ntra_bio IF ( flag_active(jn) ) THEN WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn) WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn, jk), & jk = 1, nlay_bio ) ENDIF ! flag_active END DO ! jn WRITE(numout,*) ENDIF ! ln_write_bio IF ( ln_carbon ) CALL ice_carb_chem WRITE(numout,*) WRITE(numout,*) ' End of ice_bio_diff ' WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' ! !=============================================================================! !-- End of ice_bio_diff -- END