472   !! 
473   !! blank fastsinking 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 nondiatom 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  06032017  Forgotten var to init 
574   f_omarg(:,:) = 0.0 !! 
575   f_omcal(:,:) = 0.0 
576   xFree(:,:) = 0.0 !! state variables for ironligand 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  171116  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   !! 
 238  !! 
 239  !! Initialise arrays to zero and set up arrays for diagnostics 
 240  !! 
 241  CALL bio_medusa_init( kt ) 
 242  
1724   f_wind = wndm(ji,jj) 
1725   IF (lk_oasis) THEN 
1726   f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 
1727   ENDIF 
1728   !! 
1729   !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
1730   !! in MEDUSA, the gas transfer velocity used in the carbon 
1731   !! and oxygen cycles has been harmonised and is calculated 
1732   !! by the same function here; this harmonisation includes 
1733   !! changes to the PML carbonate chemistry scheme so that 
1734   !! it too makes use of the same gas transfer velocity; the 
1735   !! preferred parameterisation of this is Wanninkhof (2014), 
1736   !! option 7 
1737   !! 
1738   # if defined key_debug_medusa 
1739   IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
1740   CALL flush(numout) 
1741   # endif 
1742   CALL gas_transfer( f_wind, 1, 7, & ! inputs 
1743   f_kw660 ) ! outputs 
1744   # if defined key_debug_medusa 
1745   IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
1746   CALL flush(numout) 
1747   # endif 
1748   !! 
1749   !! air pressure (atm); ultimately this will use air pressure at the base 
1750   !! of the UKESM1 atmosphere 
1751   !! 
1752   f_pp0 = 1.0 
1753   !! 
1754   !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp 
1755   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
1756   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 
1757   !! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind =', f_wind 
1758   !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) 
1759   !! 
1760   # if defined key_axy_carbchem 
1761   # if defined key_mocsy 
1762   !! 
1763   !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY2 carbonate 
1764   !! chemistry package; note that depth is set to 
1765   !! zero in this call 
1766   CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 
1767   f_pp0, 0.0, gphit(ji,jj), f_kw660, f_xco2a, 1, & ! inputs 
1768   f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 
1769   f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 
1770   f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 
1771   f_co2starair, f_co2flux, f_dpco2 ) ! outputs 
1772   !! 
1773   f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 > umol / kg 
1774   f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 > ueq / kg 
1775   f_dcf = f_rhosw 
1776   # else 
1777   iters = 0 
1778   !! 
1779   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP2, but not) 
1780   CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a, & ! inputs 
1781   f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 
1782   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters ) ! outputs 
1783   !! 
1784   !! AXY (09/01/14): removed iteration and NaN checks; these have 
1785   !! been moved to trc_co2_medusa together with a 
1786   !! fudge that amends erroneous values (this is 
1787   !! intended to be a temporary fudge!); the 
1788   !! output warnings are retained here so that 
1789   !! failure position can be determined 
1790   if (iters .eq. 25) then 
1791   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & 
1792   iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
1793   endif 
1794   # endif 
1795   # else 
1796   !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
1797   !! quasisensible alternatives 
1798   f_ph = 8.1 
1799   f_pco2w = f_xco2a 
1800   f_h2co3 = 0.005 * zdic 
1801   f_hco3 = 0.865 * zdic 
1802   f_co3 = 0.130 * zdic 
1803   f_omcal(ji,jj) = 4. 
1804   f_omarg(ji,jj) = 2. 
1805   f_co2flux = 0. 
1806   f_TDIC = zdic 
1807   f_TALK = zalk 
1808   f_dcf = 1.026 
1809   f_henry = 1. 
1810   !! AXY (23/06/15): add in some extra MOCSY diagnostics 
1811   f_fco2w = f_xco2a 
1812   f_BetaD = 1. 
1813   f_rhosw = 1.026 
1814   f_opres = 0. 
1815   f_insitut = ztmp 
1816   f_pco2atm = f_xco2a 
1817   f_fco2atm = f_xco2a 
1818   f_schmidtco2 = 660. 
1819   f_kwco2 = 0. 
1820   f_K0 = 0. 
1821   f_co2starair = f_xco2a 
1822   f_dpco2 = 0. 
1823   # endif 
1824   !! 
1825   !! mmol/m2/s > mmol/m3/d; correct for seaice; divide through by layer thickness 
1826   f_co2flux = (1.  fr_i(ji,jj)) * f_co2flux * 86400. / fthk 
1827   !! 
1828   !! oxygen (O2); OCMIP2 code 
1829   !! AXY (23/06/15): amend input list for oxygen to account for common gas 
1830   !! transfer velocity 
1831   !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs 
1832   !! f_kw660, f_o2flux, f_o2sat ) ! outputs 
1833   CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs 
1834   f_kwo2, f_o2flux, f_o2sat ) ! outputs 
1835   !! 
1836   !! mmol/m2/s > mol/m3/d; correct for seaice; divide through by layer thickness 
1837   f_o2flux = (1.  fr_i(ji,jj)) * f_o2flux * 86400. / fthk 
1838   !! 
1839   !! Jpalm (082014) 
1840   !! DMS surface concentration calculation 
1841   !! initialy added for UKESM1 model. 
1842   !! using METOFFICE subroutine. 
1843   !! DMS module only needs Chl concentration and MLD 
1844   !! to get an aproximate value of DMS concentration. 
1845   !! airsea fluxes are calculated by atmospheric chemitry model 
1846   !! from atm and ocsurface concentrations. 
1847   !! 
1848   !! AXY (13/03/15): this is amended to calculate all of the DMS 
1849   !! estimates examined during UKESM1 (see comments 
1850   !! in trcdms_medusa.F90) 
1851   !! 
1852   !! AXY (25/05/17): amended to additionally pass DIN limitation as well as [DIN]; 
1853   !! accounts for differences in nutrient halfsaturations; changes 
1854   !! also made in trc_dms_medusa; this permits an additional DMS 
1855   !! calculation while retaining the existing Anderson one 
1856   !! 
1857   IF (jdms .eq. 1) THEN 
1858   !! 
1859   !! calculate weighted halfsaturation for DIN uptake 
1860   dms_wtkn = ((zphn * xnln) + (zphd * xnld)) / (zphn + zphd) 
1861   !! 
1862   !! feed in correct inputs 
1863   if (jdms_input .eq. 0) then 
1864   !! use instantaneous inputs 
1865   dms_nlim = zdin / (zdin + dms_wtkn) 
1866   !! 
1867   CALL trc_dms_medusa( zchn, zchd, & ! inputs 
1868   hmld(ji,jj), qsr(ji,jj), & ! inputs 
1869   zdin, dms_nlim, & ! inputs 
1870   dms_andr, dms_simo, dms_aran, dms_hall, dms_andm ) ! outputs 
1871   else 
1872   !! use dielaverage inputs 
1873   dms_nlim = zn_dms_din(ji,jj) / (zn_dms_din(ji,jj) + dms_wtkn) 
1874   !! 
1875   CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs 
1876   zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), & ! inputs 
1877   zn_dms_din(ji,jj), dms_nlim, & ! inputs 
1878   dms_andr, dms_simo, dms_aran, dms_hall, dms_andm ) ! outputs 
1879   endif 
1880   !! 
1881   !! assign correct output to variable passed to atmosphere 
1882   if (jdms_model .eq. 1) then 
1883   dms_surf = dms_andr 
1884   elseif (jdms_model .eq. 2) then 
1885   dms_surf = dms_simo 
1886   elseif (jdms_model .eq. 3) then 
1887   dms_surf = dms_aran 
1888   elseif (jdms_model .eq. 4) then 
1889   dms_surf = dms_hall 
1890   elseif (jdms_model .eq. 5) then 
1891   dms_surf = dms_andm 
1892   endif 
1893   !! 
1894   !! 2D diag through iom_use 
1895   IF( lk_iomput ) THEN 
1896   IF( med_diag%DMS_SURF%dgsave ) THEN 
1897   dms_surf2d(ji,jj) = dms_surf 
1898   ENDIF 
1899   IF( med_diag%DMS_ANDR%dgsave ) THEN 
1900   dms_andr2d(ji,jj) = dms_andr 
1901   ENDIF 
1902   IF( med_diag%DMS_SIMO%dgsave ) THEN 
1903   dms_simo2d(ji,jj) = dms_simo 
1904   ENDIF 
1905   IF( med_diag%DMS_ARAN%dgsave ) THEN 
1906   dms_aran2d(ji,jj) = dms_aran 
1907   ENDIF 
1908   IF( med_diag%DMS_HALL%dgsave ) THEN 
1909   dms_hall2d(ji,jj) = dms_hall 
1910   ENDIF 
1911   IF( med_diag%DMS_ANDM%dgsave ) THEN 
1912   dms_andm2d(ji,jj) = dms_andm 
1913   ENDIF 
1914   # if defined key_debug_medusa 
1915   IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' 
1916   CALL flush(numout) 
1917   # endif 
1918   ENDIF 
1919   !! End iom 
1920   ENDIF 
1921   !! End DMS Loop 
1922   !! 
1923   !! store 2D outputs 
1924   !! 
1925   !! JPALM  171116  put fgco2 out of diag request 
1926   !! is needed for coupling; pass through restart 
1927   !! IF( med_diag%FGCO2%dgsave ) THEN 
1928   !! convert from mol/m2/day to kg/m2/s 
1929   fgco2(ji,jj) = f_co2flux * fthk * CO2flux_conv !! mmolC/m3/d > kgCO2/m2/s 
1930   !! ENDIF 
1931   IF ( lk_iomput ) THEN 
1932   IF( med_diag%ATM_PCO2%dgsave ) THEN 
1933   f_pco2a2d(ji,jj) = f_pco2atm 
1934   ENDIF 
1935   IF( med_diag%OCN_PCO2%dgsave ) THEN 
1936   f_pco2w2d(ji,jj) = f_pco2w 
1937   ENDIF 
1938   IF( med_diag%CO2FLUX%dgsave ) THEN 
1939   f_co2flux2d(ji,jj) = f_co2flux * fthk !! mmol/m3/d > mmol/m2/d 
1940   ENDIF 
1941   IF( med_diag%TCO2%dgsave ) THEN 
1942   f_TDIC2d(ji,jj) = f_TDIC 
1943   ENDIF 
1944   IF( med_diag%TALK%dgsave ) THEN 
1945   f_TALK2d(ji,jj) = f_TALK 
1946   ENDIF 
1947   IF( med_diag%KW660%dgsave ) THEN 
1948   f_kw6602d(ji,jj) = f_kw660 
1949   ENDIF 
1950   IF( med_diag%ATM_PP0%dgsave ) THEN 
1951   f_pp02d(ji,jj) = f_pp0 
1952   ENDIF 
1953   IF( med_diag%O2FLUX%dgsave ) THEN 
1954   f_o2flux2d(ji,jj) = f_o2flux 
1955   ENDIF 
1956   IF( med_diag%O2SAT%dgsave ) THEN 
1957   f_o2sat2d(ji,jj) = f_o2sat 
1958   ENDIF 
1959   !! AXY (24/11/16): add in extra MOCSY diagnostics 
1960   IF( med_diag%ATM_XCO2%dgsave ) THEN 
1961   f_xco2a_2d(ji,jj) = f_xco2a 
1962   ENDIF 
1963   IF( med_diag%OCN_FCO2%dgsave ) THEN 
1964   f_fco2w_2d(ji,jj) = f_fco2w 
1965   ENDIF 
1966   IF( med_diag%ATM_FCO2%dgsave ) THEN 
1967   f_fco2a_2d(ji,jj) = f_fco2atm 
1968   ENDIF 
1969   IF( med_diag%OCN_RHOSW%dgsave ) THEN 
1970   f_ocnrhosw_2d(ji,jj) = f_rhosw 
1971   ENDIF 
1972   IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
1973   f_ocnschco2_2d(ji,jj) = f_schmidtco2 
1974   ENDIF 
1975   IF( med_diag%OCN_KWCO2%dgsave ) THEN 
1976   f_ocnkwco2_2d(ji,jj) = f_kwco2 
1977   ENDIF 
1978   IF( med_diag%OCN_K0%dgsave ) THEN 
1979   f_ocnk0_2d(ji,jj) = f_K0 
1980   ENDIF 
1981   IF( med_diag%CO2STARAIR%dgsave ) THEN 
1982   f_co2starair_2d(ji,jj) = f_co2starair 
1983   ENDIF 
1984   IF( med_diag%OCN_DPCO2%dgsave ) THEN 
1985   f_ocndpco2_2d(ji,jj) = f_dpco2 
1986   ENDIF 
1987   ENDIF 
1988   !! 
1989   endif 
1990   !! End jk = 1 loop within ROAM key 
1991   
1992   !! AXY (11/11/16): CMIP6 oxygen saturation 3D diagnostic 
1993   IF ( med_diag%O2SAT3%dgsave ) THEN 
1994   call oxy_sato( ztmp, zsal, f_o2sat3 ) 
1995   o2sat3(ji, jj, jk) = f_o2sat3 
1996   ENDIF 
1997   
1998   # endif 
1999   
2000   if ( jk .eq. 1 ) then 
2001   !! 
2002   !! River inputs 
2003   !! 
2004   !! 
2005   !! runoff comes in as kg / m2 / s 
2006   !! used and written out as m3 / m2 / d (= m / d) 
2007   !! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 
2008   !! 
2009   !! AXY (17/07/14): the compiler doesn't like this line for some reason; 
2010   !! as MEDUSA doesn't even use runoff for riverine inputs, 
2011   !! a temporary solution is to switch off runoff entirely 
2012   !! here; again, this change is one of several that will 
2013   !! need revisiting once MEDUSA has bedded down in UKESM1; 
2014   !! particularly so if the land scheme provides information 
2015   !! concerning nutrient fluxes 
2016   !! 
2017   !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 
2018   f_runoff(ji,jj) = 0.0 
2019   !! 
2020   !! nutrients are added via rivers to the model in one of two ways: 
2021   !! 1. via river concentration; i.e. the average nutrient concentration 
2022   !! of a river water is described by a spatial file, and this is 
2023   !! multiplied by runoff to give a nutrient flux 
2024   !! 2. via direct river flux; i.e. the average nutrient flux due to 
2025   !! rivers is described by a spatial file, and this is simply applied 
2026   !! as a direct nutrient flux (i.e. it does not relate or respond to 
2027   !! model runoff) 
2028   !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 
2029   !! alkalinity are derived from continentscale DIC estimates (Huang et al., 
2030   !! 2012) and some Arctic river alkalinity estimates (Katya?) 
2031   !! 
2032   !! as of 19/07/12, riverine nutrients can now be spread vertically across 
2033   !! several grid cells rather than just poured into the surface box; this 
2034   !! block of code is still executed, however, to set up the total amounts 
2035   !! of nutrient entering via rivers 
2036   !! 
2037   !! nitrogen 
2038   if (jriver_n .eq. 1) then 
2039   !! river concentration specified; use runoff to calculate input 
2040   f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 
2041   elseif (jriver_n .eq. 2) then 
2042   !! river flux specified; independent of runoff 
2043   f_riv_n(ji,jj) = riv_n(ji,jj) 
2044   endif 
2045   !! 
2046   !! silicon 
2047   if (jriver_si .eq. 1) then 
2048   !! river concentration specified; use runoff to calculate input 
2049   f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 
2050   elseif (jriver_si .eq. 2) then 
2051   !! river flux specified; independent of runoff 
2052   f_riv_si(ji,jj) = riv_si(ji,jj) 
2053   endif 
2054   !! 
2055   !! carbon 
2056   if (jriver_c .eq. 1) then 
2057   !! river concentration specified; use runoff to calculate input 
2058   f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 
2059   elseif (jriver_c .eq. 2) then 
2060   !! river flux specified; independent of runoff 
2061   f_riv_c(ji,jj) = riv_c(ji,jj) 
2062   endif 
2063   !! 
2064   !! alkalinity 
2065   if (jriver_alk .eq. 1) then 
2066   !! river concentration specified; use runoff to calculate input 
2067   f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 
2068   elseif (jriver_alk .eq. 2) then 
2069   !! river flux specified; independent of runoff 
2070   f_riv_alk(ji,jj) = riv_alk(ji,jj) 
2071   endif 
2072   
2073   endif 
2074   
2075   !! 
2076   !! Chlorophyll calculations 
2077   !! 
2078   !! 
2079   !! nondiatoms 
2080   if (zphn.GT.rsmall) then 
2081   fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn))) 
2082   faln = xaln * fthetan 
2083   else 
2084   fthetan = 0. 
2085   faln = 0. 
2086   endif 
2087   !! 
2088   !! diatoms 
2089   if (zphd.GT.rsmall) then 
2090   fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd))) 
2091   fald = xald * fthetad 
2092   else 
2093   fthetad = 0. 
2094   fald = 0. 
2095   endif 
2096   
2097   # if defined key_debug_medusa 
2098   !! report biological calculations 
2099   if (idf.eq.1.AND.idfval.eq.1) then 
2100   IF (lwp) write (numout,*) '' 
2101   IF (lwp) write (numout,*) 'faln(',jk,') = ', faln 
2102   IF (lwp) write (numout,*) 'fald(',jk,') = ', fald 
2103   endif 
2104   # endif 
2105   
2106   !! 
2107   !! Phytoplankton light limitation 
2108   !! 
2109   !! 
2110   !! It is assumed xpar is the depthaveraged (vertical layer) PAR 
2111   !! Light limitation (check selfshading) in W/m2 
2112   !! 
2113   !! Note that there is no temperature dependence in phytoplankton 
2114   !! growth rate or any other function. 
2115   !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid 
2116   !! NaNs in case of Phy==0. 
2117   !! 
2118   !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and nondiat: 
2119   !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012 
2120   !! 
2121   !! AXY (16/07/09) 
2122   !! temperature for new Eppley style phytoplankton growth 
2123   loc_T = tsn(ji,jj,jk,jp_tem) 
2124   fun_T = 1.066**(1.0 * loc_T) 
2125   !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for 
2126   !phytoplankton 
2127   !! growth; remin. unaffected 
2128   fun_Q10 = jq10**((loc_T  0.0) / 10.0) 
2129   if (jphy.eq.1) then 
2130   xvpnT = xvpn * fun_T 
2131   xvpdT = xvpd * fun_T 
2132   elseif (jphy.eq.2) then 
2133   xvpnT = xvpn * fun_Q10 
2134   xvpdT = xvpd * fun_Q10 
2135   else 
2136   xvpnT = xvpn 
2137   xvpdT = xvpd 
2138   endif 
2139   !! 
2140   !! nondiatoms 
2141   fchn1 = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
2142   if (fchn1.GT.rsmall) then 
2143   fchn = xvpnT / (sqrt(fchn1) + tiny(fchn1)) 
2144   else 
2145   fchn = 0. 
2146   endif 
2147   fjln = fchn * faln * xpar(ji,jj,jk) !! nondiatom J term 
2148   fjlim_pn = fjln / xvpnT 
2149   !! 
2150   !! diatoms 
2151   fchd1 = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
2152   if (fchd1.GT.rsmall) then 
2153   fchd = xvpdT / (sqrt(fchd1) + tiny(fchd1)) 
2154   else 
2155   fchd = 0. 
2156   endif 
2157   fjld = fchd * fald * xpar(ji,jj,jk) !! diatom J term 
2158   fjlim_pd = fjld / xvpdT 
2159   
2160   # if defined key_debug_medusa 
2161   !! report phytoplankton light limitation 
2162   if (idf.eq.1.AND.idfval.eq.1) then 
2163   IF (lwp) write (numout,*) '' 
2164   IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn 
2165   IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd 
2166   IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln 
2167   IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld 
2168   endif 
2169   # endif 
2170   
2171   !! 
2172   !! Phytoplankton nutrient limitation 
2173   !! 
2174   !! 
2175   !! nondiatoms (N, Fe) 
2176   fnln = zdin / (zdin + xnln) !! nondiatom Qn term 
2177   ffln = zfer / (zfer + xfln) !! nondiatom Qf term 
2178   !! 
2179   !! diatoms (N, Si, Fe) 
2180   fnld = zdin / (zdin + xnld) !! diatom Qn term 
2181   fsld = zsil / (zsil + xsld) !! diatom Qs term 
2182   ffld = zfer / (zfer + xfld) !! diatom Qf term 
2183   
2184   # if defined key_debug_medusa 
2185   !! report phytoplankton nutrient limitation 
2186   if (idf.eq.1.AND.idfval.eq.1) then 
2187   IF (lwp) write (numout,*) '' 
2188   IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln 
2189   IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld 
2190   IF (lwp) write (numout,*) 'ffln(',jk,') = ', ffln 
2191   IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld 
2192   IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld 
2193   endif 
2194   # endif 
2195   
2196   !! 
2197   !! Primary production (nondiatoms) 
2198   !! (note: still needs multiplying by phytoplankton concentration) 
2199   !! 
2200   !! 
2201   if (jliebig .eq. 0) then 
2202   !! multiplicative nutrient limitation 
2203   fpnlim = fnln * ffln 
2204   elseif (jliebig .eq. 1) then 
2205   !! Liebig Law (= most limiting) nutrient limitation 
2206   fpnlim = min(fnln, ffln) 
2207   endif 
2208   fprn = fjln * fpnlim 
2209   
2210   !! 
2211   !! Primary production (diatoms) 
2212   !! (note: still needs multiplying by phytoplankton concentration) 
2213   !! 
2214   !! production here is split between nitrogen production and that of 
2215   !! silicon; depending upon the "intracellular" ratio of Si:N, model 
2216   !! diatoms will uptake nitrogen/silicon differentially; this borrows 
2217   !! from the diatom model of Mongin et al. (2006) 
2218   !! 
2219   !! 
2220   if (jliebig .eq. 0) then 
2221   !! multiplicative nutrient limitation 
2222   fpdlim = fnld * ffld 
2223   elseif (jliebig .eq. 1) then 
2224   !! Liebig Law (= most limiting) nutrient limitation 
2225   fpdlim = min(fnld, ffld) 
2226   endif 
2227   !! 
2228   if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then 
2229   !! "intracellular" elemental ratios 
2230   ! fsin = zpds / (zphd + tiny(zphd)) 
2231   ! fnsi = zphd / (zpds + tiny(zpds)) 
2232   fsin = 0.0 
2233   IF( zphd .GT. rsmall) fsin = zpds / zphd 
2234   fnsi = 0.0 
2235   IF( zpds .GT. rsmall) fnsi = zphd / zpds 
2236   !! AXY (23/02/10): these next variables derive from Mongin et al. (2003) 
2237   fsin1 = 3.0 * xsin0 !! = 0.6 
2238   fnsi1 = 1.0 / fsin1 !! = 1.667 
2239   fnsi2 = 1.0 / xsin0 !! = 5.0 
2240   !! 
2241   !! conditionalities based on ratios 
2242   !! nitrogen (and iron and carbon) 
2243   if (fsin.le.xsin0) then 
2244   fprd = 0.0 
2245   fsld2 = 0.0 
2246   elseif (fsin.lt.fsin1) then 
2247   fprd = xuif * ((fsin  xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim) 
2248   fsld2 = xuif * ((fsin  xsin0) / (fsin + tiny(fsin))) 
2249   elseif (fsin.ge.fsin1) then 
2250   fprd = (fjld * fpdlim) 
2251   fsld2 = 1.0 
2252   endif 
2253   !! 
2254   !! silicon 
2255   if (fsin.lt.fnsi1) then 
2256   fprds = (fjld * fsld) 
2257   elseif (fsin.lt.fnsi2) then 
2258   fprds = xuif * ((fnsi  xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld) 
2259   else 
2260   fprds = 0.0 
2261   endif 
2262   else 
2263   fsin = 0.0 
2264   fnsi = 0.0 
2265   fprd = 0.0 
2266   fsld2 = 0.0 
2267   fprds = 0.0 
2268   endif 
2269   
2270   # if defined key_debug_medusa 
2271   !! report phytoplankton growth (including diatom silicon submodel) 
2272   if (idf.eq.1.AND.idfval.eq.1) then 
2273   IF (lwp) write (numout,*) '' 
2274   IF (lwp) write (numout,*) 'fsin(',jk,') = ', fsin 
2275   IF (lwp) write (numout,*) 'fnsi(',jk,') = ', fnsi 
2276   IF (lwp) write (numout,*) 'fsld2(',jk,') = ', fsld2 
2277   IF (lwp) write (numout,*) 'fprn(',jk,') = ', fprn 
2278   IF (lwp) write (numout,*) 'fprd(',jk,') = ', fprd 
2279   IF (lwp) write (numout,*) 'fprds(',jk,') = ', fprds 
2280   endif 
2281   # endif 
2282   
2283   !! 
2284   !! Mixed layer primary production 
2285   !! this block calculates the amount of primary production that occurs 
2286   !! within the upper mixed layer; this allows the separate diagnosis 
2287   !! of "subsurface" primary production; it does assume that short 
2288   !! term variability in mixed layer depth doesn't mess with things 
2289   !! though 
2290   !! 
2291   !! 
2292   if (fdep1.le.hmld(ji,jj)) then 
2293   !! this level is entirely in the mixed layer 
2294   fq0 = 1.0 
2295   elseif (fdep.ge.hmld(ji,jj)) then 
2296   !! this level is entirely below the mixed layer 
2297   fq0 = 0.0 
2298   else 
2299   !! this level straddles the mixed layer 
2300   fq0 = (hmld(ji,jj)  fdep) / fthk 
2301   endif 
2302   !! 
2303   fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0) 
2304   fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0) 
2305   
2306   !! 
2307   !! Vertical Integral  
2308   !! 
2309   ftot_pn(ji,jj) = ftot_pn(ji,jj) + (zphn * fthk) !! vertical integral nondiatom phytoplankton 
2310   ftot_pd(ji,jj) = ftot_pd(ji,jj) + (zphd * fthk) !! vertical integral diatom phytoplankton 
2311   ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk) !! vertical integral microzooplankton 
2312   ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk) !! vertical integral mesozooplankton 
2313   ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk) !! vertical integral slow detritus, nitrogen 
2314   ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk) !! vertical integral slow detritus, carbon 
2315   
2316   !! 
2317   !! More chlorophyll calculations 
2318   !! 
2319   !! 
2320   !! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk))) 
2321   !! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk))) 
2322   frn = (xthetam * fchn * fnln * ffln ) / (fthetan + tiny(fthetan)) 
2323   !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid 
2324   !! limitation, is used in the following line to regulate chlorophyll 
2325   !! growth in a manner that is inconsistent with its use in the regulation 
2326   !! of biomass growth; the Mongin term term used in growth is more complex 
2327   !! than the simple multiplicative function used below 
2328   !! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad)) 
2329   !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to 
2330   !! regulate chlorophyll growth 
2331   frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad)) 
2332   
2333   # if defined key_debug_medusa 
2334   !! report chlorophyll calculations 
2335   if (idf.eq.1.AND.idfval.eq.1) then 
2336   IF (lwp) write (numout,*) '' 
2337   IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan 
2338   IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad 
2339   IF (lwp) write (numout,*) 'frn(',jk,') = ', frn 
2340   IF (lwp) write (numout,*) 'frd(',jk,') = ', frd 
2341   endif 
2342   # endif 
2343   
2344   !! 
2345   !! Zooplankton Grazing 
2346   !! this code supplements the base grazing model with one that 
2347   !! considers the C:N ratio of grazed food and balances this against 
2348   !! the requirements of zooplankton growth; this model is derived 
2349   !! from that of Anderson & Pondaven (2003) 
2350   !! 
2351   !! the current version of the code assumes a fixed C:N ratio for 
2352   !! detritus (in contrast to Anderson & Pondaven, 2003), though the 
2353   !! full equations are retained for future extension 
2354   !! 
2355   !! 
2356   !! 
2357   !! Microzooplankton first 
2358   !! 
2359   !! 
2360   fmi1 = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet) 
2361   fmi = xgmi * zzmi / fmi1 
2362   fgmipn = fmi * xpmipn * zphn * zphn !! grazing on nondiatoms 
2363   fgmid = fmi * xpmid * zdet * zdet !! grazing on detrital nitrogen 
2364   # if defined key_roam 
2365   fgmidc = rsmall !acc 
2366   IF ( zdet .GT. rsmall ) fgmidc = (zdtc / (zdet + tiny(zdet))) * fgmid !! grazing on detrital carbon 
2367   # else 
2368   !! AXY (26/11/08): implicit detrital carbon change 
2369   fgmidc = xthetad * fgmid !! grazing on detrital carbon 
2370   # endif 
2371   !! 
2372   !! which translates to these incoming N and C fluxes 
2373   finmi = (1.0  xphi) * (fgmipn + fgmid) 
2374   ficmi = (1.0  xphi) * ((xthetapn * fgmipn) + fgmidc) 
2375   !! 
2376   !! the ideal food C:N ratio for microzooplankton 
2377   !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
2378   fstarmi = (xbetan * xthetazmi) / (xbetac * xkc) 
2379   !! 
2380   !! process these to determine proportioning of grazed N and C 
2381   !! (since there is no explicit consideration of respiration, 
2382   !! only growth and excretion are calculated here) 
2383   fmith = (ficmi / (finmi + tiny(finmi))) 
2384   if (fmith.ge.fstarmi) then 
2385   fmigrow = xbetan * finmi 
2386   fmiexcr = 0.0 
2387   else 
2388   fmigrow = (xbetac * xkc * ficmi) / xthetazmi 
2389   fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith)))  ((xbetac * xkc) / xthetazmi)) 
2390   endif 
2391   # if defined key_roam 
2392   fmiresp = (xbetac * ficmi)  (xthetazmi * fmigrow) 
2393   # endif 
2394   
2395   # if defined key_debug_medusa 
2396   !! report microzooplankton grazing 
2397   if (idf.eq.1.AND.idfval.eq.1) then 
2398   IF (lwp) write (numout,*) '' 
2399   IF (lwp) write (numout,*) 'fmi1(',jk,') = ', fmi1 
2400   IF (lwp) write (numout,*) 'fmi(',jk,') = ', fmi 
2401   IF (lwp) write (numout,*) 'fgmipn(',jk,') = ', fgmipn 
2402   IF (lwp) write (numout,*) 'fgmid(',jk,') = ', fgmid 
2403   IF (lwp) write (numout,*) 'fgmidc(',jk,') = ', fgmidc 
2404   IF (lwp) write (numout,*) 'finmi(',jk,') = ', finmi 
2405   IF (lwp) write (numout,*) 'ficmi(',jk,') = ', ficmi 
2406   IF (lwp) write (numout,*) 'fstarmi(',jk,') = ', fstarmi 
2407   IF (lwp) write (numout,*) 'fmith(',jk,') = ', fmith 
2408   IF (lwp) write (numout,*) 'fmigrow(',jk,') = ', fmigrow 
2409   IF (lwp) write (numout,*) 'fmiexcr(',jk,') = ', fmiexcr 
2410   # if defined key_roam 
2411   IF (lwp) write (numout,*) 'fmiresp(',jk,') = ', fmiresp 
2412   # endif 
2413   endif 
2414   # endif 
2415   
2416   !! 
2417   !! Mesozooplankton second 
2418   !! 
2419   !! 
2420   fme1 = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + & 
2421   (xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet) 
2422   fme = xgme * zzme / fme1 
2423   fgmepn = fme * xpmepn * zphn * zphn !! grazing on nondiatoms 
2424   fgmepd = fme * xpmepd * zphd * zphd !! grazing on diatoms 
2425   fgmepds = fsin * fgmepd !! grazing on diatom silicon 
2426   fgmezmi = fme * xpmezmi * zzmi * zzmi !! grazing on microzooplankton 
2427   fgmed = fme * xpmed * zdet * zdet !! grazing on detrital nitrogen 
2428   # if defined key_roam 
2429   fgmedc = rsmall !acc 
2430   IF ( zdet .GT. rsmall ) fgmedc = (zdtc / (zdet + tiny(zdet))) * fgmed !! grazing on detrital carbon 
2431   # else 
2432   !! AXY (26/11/08): implicit detrital carbon change 
2433   fgmedc = xthetad * fgmed !! grazing on detrital carbon 
2434   # endif 
2435   !! 
2436   !! which translates to these incoming N and C fluxes 
2437   finme = (1.0  xphi) * (fgmepn + fgmepd + fgmezmi + fgmed) 
2438   ficme = (1.0  xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & 
2439   (xthetazmi * fgmezmi) + fgmedc) 
2440   !! 
2441   !! the ideal food C:N ratio for mesozooplankton 
2442   !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
2443   fstarme = (xbetan * xthetazme) / (xbetac * xkc) 
2444   !! 
2445   !! process these to determine proportioning of grazed N and C 
2446   !! (since there is no explicit consideration of respiration, 
2447   !! only growth and excretion are calculated here) 
2448   fmeth = (ficme / (finme + tiny(finme))) 
2449   if (fmeth.ge.fstarme) then 
2450   fmegrow = xbetan * finme 
2451   fmeexcr = 0.0 
2452   else 
2453   fmegrow = (xbetac * xkc * ficme) / xthetazme 
2454   fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth)))  ((xbetac * xkc) / xthetazme)) 
2455   endif 
2456   # if defined key_roam 
2457   fmeresp = (xbetac * ficme)  (xthetazme * fmegrow) 
2458   # endif 
2459   
2460   # if defined key_debug_medusa 
2461   !! report mesozooplankton grazing 
2462   if (idf.eq.1.AND.idfval.eq.1) then 
2463   IF (lwp) write (numout,*) '' 
2464   IF (lwp) write (numout,*) 'fme1(',jk,') = ', fme1 
2465   IF (lwp) write (numout,*) 'fme(',jk,') = ', fme 
2466   IF (lwp) write (numout,*) 'fgmepn(',jk,') = ', fgmepn 
2467   IF (lwp) write (numout,*) 'fgmepd(',jk,') = ', fgmepd 
2468   IF (lwp) write (numout,*) 'fgmepds(',jk,') = ', fgmepds 
2469   IF (lwp) write (numout,*) 'fgmezmi(',jk,') = ', fgmezmi 
2470   IF (lwp) write (numout,*) 'fgmed(',jk,') = ', fgmed 
2471   IF (lwp) write (numout,*) 'fgmedc(',jk,') = ', fgmedc 
2472   IF (lwp) write (numout,*) 'finme(',jk,') = ', finme 
2473   IF (lwp) write (numout,*) 'ficme(',jk,') = ', ficme 
2474   IF (lwp) write (numout,*) 'fstarme(',jk,') = ', fstarme 
2475   IF (lwp) write (numout,*) 'fmeth(',jk,') = ', fmeth 
2476   IF (lwp) write (numout,*) 'fmegrow(',jk,') = ', fmegrow 
2477   IF (lwp) write (numout,*) 'fmeexcr(',jk,') = ', fmeexcr 
2478   # if defined key_roam 
2479   IF (lwp) write (numout,*) 'fmeresp(',jk,') = ', fmeresp 
2480   # endif 
2481   endif 
2482   # endif 
2483   
2484   fzmi_i(ji,jj) = fzmi_i(ji,jj) + fthk * ( & 
2485   fgmipn + fgmid ) 
2486   fzmi_o(ji,jj) = fzmi_o(ji,jj) + fthk * ( & 
2487   fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0  xbetan) * finmi) ) 
2488   fzme_i(ji,jj) = fzme_i(ji,jj) + fthk * ( & 
2489   fgmepn + fgmepd + fgmezmi + fgmed ) 
2490   fzme_o(ji,jj) = fzme_o(ji,jj) + fthk * ( & 
2491   fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0  xbetan) * finme) ) 
2492   
2493   !! 
2494   !! Plankton metabolic losses 
2495   !! Linear loss processes assumed to be metabolic in origin 
2496   !! 
2497   !! 
2498   fdpn2 = xmetapn * zphn 
2499   fdpd2 = xmetapd * zphd 
2500   fdpds2 = xmetapd * zpds 
2501   fdzmi2 = xmetazmi * zzmi 
2502   fdzme2 = xmetazme * zzme 
2503   
2504   !! 
2505   !! Plankton mortality losses 
2506   !! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced 
2507   !! to improve performance in gyres 
2508   !! 
2509   !! 
2510   !! nondiatom phytoplankton 
2511   if (jmpn.eq.1) fdpn = xmpn * zphn !! linear 
2512   if (jmpn.eq.2) fdpn = xmpn * zphn * zphn !! quadratic 
2513   if (jmpn.eq.3) fdpn = xmpn * zphn * & !! hyperbolic 
2514   (zphn / (xkphn + zphn)) 
2515   if (jmpn.eq.4) fdpn = xmpn * zphn * & !! sigmoid 
2516   ((zphn * zphn) / (xkphn + (zphn * zphn))) 
2517   !! 
2518   !! diatom phytoplankton 
2519   if (jmpd.eq.1) fdpd = xmpd * zphd !! linear 
2520   if (jmpd.eq.2) fdpd = xmpd * zphd * zphd !! quadratic 
2521   if (jmpd.eq.3) fdpd = xmpd * zphd * & !! hyperbolic 
2522   (zphd / (xkphd + zphd)) 
2523   if (jmpd.eq.4) fdpd = xmpd * zphd * & !! sigmoid 
2524   ((zphd * zphd) / (xkphd + (zphd * zphd))) 
2525   fdpds = fdpd * fsin 
2526   !! 
2527   !! microzooplankton 
2528   if (jmzmi.eq.1) fdzmi = xmzmi * zzmi !! linear 
2529   if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi !! quadratic 
2530   if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * & !! hyperbolic 
2531   (zzmi / (xkzmi + zzmi)) 
2532   if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * & !! sigmoid 
2533   ((zzmi * zzmi) / (xkzmi + (zzmi * zzmi))) 
2534   !! 
2535   !! mesozooplankton 
2536   if (jmzme.eq.1) fdzme = xmzme * zzme !! linear 
2537   if (jmzme.eq.2) fdzme = xmzme * zzme * zzme !! quadratic 
2538   if (jmzme.eq.3) fdzme = xmzme * zzme * & !! hyperbolic 
2539   (zzme / (xkzme + zzme)) 
2540   if (jmzme.eq.4) fdzme = xmzme * zzme * & !! sigmoid 
2541   ((zzme * zzme) / (xkzme + (zzme * zzme))) 
2542   
2543   !! 
2544   !! Detritus remineralisation 
2545   !! Constant or temperaturedependent 
2546   !! 
2547   !! 
2548   if (jmd.eq.1) then 
2549   !! temperaturedependent 
2550   fdd = xmd * fun_T * zdet 
2551   # if defined key_roam 
2552   fddc = xmdc * fun_T * zdtc 
2553   # endif 
2554   elseif (jmd.eq.2) then 
2555   !! AXY (16/05/13): add in Q10based parameterisation (def in nmlst) 
2556   !! temperaturedependent 
2557   fdd = xmd * fun_Q10 * zdet 
2558   # if defined key_roam 
2559   fddc = xmdc * fun_Q10 * zdtc 
2560   # endif 
2561   else 
2562   !! temperatureindependent 
2563   fdd = xmd * zdet 
2564   # if defined key_roam 
2565   fddc = xmdc * zdtc 
2566   # endif 
2567   endif 
2568   !! 
2569   !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box 
2570   if ((jk.eq.jmbathy) .and. jsfd.eq.1) then 
2571   fdd = 1.0 * zdet 
2572   # if defined key_roam 
2573   fddc = 1.0 * zdtc 
2574   # endif 
2575   endif 
2576   
2577   # if defined key_debug_medusa 
2578   !! report plankton mortality and remineralisation 
2579   if (idf.eq.1.AND.idfval.eq.1) then 
2580   IF (lwp) write (numout,*) '' 
2581   IF (lwp) write (numout,*) 'fdpn2(',jk,') = ', fdpn2 
2582   IF (lwp) write (numout,*) 'fdpd2(',jk,') = ', fdpd2 
2583   IF (lwp) write (numout,*) 'fdpds2(',jk,')= ', fdpds2 
2584   IF (lwp) write (numout,*) 'fdzmi2(',jk,')= ', fdzmi2 
2585   IF (lwp) write (numout,*) 'fdzme2(',jk,')= ', fdzme2 
2586   IF (lwp) write (numout,*) 'fdpn(',jk,') = ', fdpn 
2587   IF (lwp) write (numout,*) 'fdpd(',jk,') = ', fdpd 
2588   IF (lwp) write (numout,*) 'fdpds(',jk,') = ', fdpds 
2589   IF (lwp) write (numout,*) 'fdzmi(',jk,') = ', fdzmi 
2590   IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme 
2591   IF (lwp) write (numout,*) 'fdd(',jk,') = ', fdd 
2592   # if defined key_roam 
2593   IF (lwp) write (numout,*) 'fddc(',jk,') = ', fddc 
2594   # endif 
2595   endif 
2596   # endif 
2597   
2598   !! 
2599   !! Detritus addition to benthos 
2600   !! If activated, slow detritus in the bottom box will enter the 
2601   !! benthic pool 
2602   !! 
2603   !! 
2604   if ((jk.eq.jmbathy) .and. jorgben.eq.1) then 
2605   !! this is the BOTTOM OCEAN BOX > into the benthic pool! 
2606   !! 
2607   f_sbenin_n(ji,jj) = (zdet * vsed * 86400.) 
2608   f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn) 
2609   # if defined key_roam 
2610   f_sbenin_c(ji,jj) = (zdtc * vsed * 86400.) 
2611   # else 
2612   f_sbenin_c(ji,jj) = (zdet * vsed * 86400. * xthetad) 
2613   # endif 
2614   endif 
2615   
2616   !! 
2617   !! Iron chemistry and fractionation 
2618   !! following the Parekh et al. (2004) scheme adopted by the Met. 
2619   !! Office, Medusa models total iron but considers "free" and 
2620   !! ligandbound forms for the purposes of scavenging (only "free" 
2621   !! iron can be scavenged 
2622   !! 
2623   !! 
2624   !! total iron concentration (mmol Fe / m3 > umol Fe / m3) 
2625   xFeT = zfer * 1.e3 
2626   !! 
2627   !! calculate fractionation (based on DiatHadOCC; in turn based on Parekh et al., 2004) 
2628   xb_coef_tmp = xk_FeL * (xLgT  xFeT)  1.0 
2629   xb2M4ac = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0) 
2630   !! 
2631   !! "free" ligand concentration 
2632   xLgF = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL 
2633   !! 
2634   !! ligandbound iron concentration 
2635   xFeL = xLgT  xLgF 
2636   !! 
2637   !! "free" iron concentration (and convert to mmol Fe / m3) 
2638   xFeF = (xFeT  xFeL) * 1.e3 
2639   xFree(ji,jj)= xFeF / (zfer + tiny(zfer)) 
2640   !! 
2641   !! scavenging of iron (multiple schemes); I'm only really happy with the 
2642   !! first one at the moment  the others involve assumptions (sometimes 
2643   !! guessed at by me) that are potentially questionable 
2644   !! 
2645   if (jiron.eq.1) then 
2646   !! 
2647   !! Scheme 1: Dutkiewicz et al. (2005) 
2648   !! This scheme includes a single scavenging term based solely on a 
2649   !! fixed rate and the availablility of "free" iron 
2650   !! 
2651   !! 
2652   ffescav = xk_sc_Fe * xFeF ! = mmol/m3/d 
2653   !! 
2654   !! 
2655   !! 
2656   !! Mick's code contains a further (optional) implicit "scavenging" of 
2657   !! iron that sets an upper bound on "free" iron concentration, and 
2658   !! essentially caps the concentration of total iron as xFeL + "free" 
2659   !! iron; since the former is constrained by a fixed total ligand 
2660   !! concentration (= 1.0 umol/m3), and the latter isn't allowed above 
2661   !! this upper bound, total iron is constrained to a maximum of ... 
2662   !! 
2663   !! xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3 
2664   !! 
2665   !! In Mick's code, the actual value of total iron is reset to this 
2666   !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't 
2667   !! our favoured approach to tracer updating here (not least because 
2668   !! of the leapfrog), so here the amount scavenged is augmented by an 
2669   !! additional amount that serves to drag total iron back towards that 
2670   !! expected from this limitation on iron concentration ... 
2671   !! 
2672   xmaxFeF = min((xFeF * 1.e3), 0.3) ! = umol/m3 
2673   !! 
2674   !! Here, the difference between current total Fe and (FeL + Fe') is 
2675   !! calculated and added to the scavenging flux already calculated 
2676   !! above ... 
2677   !! 
2678   fdeltaFe = (xFeT  (xFeL + xmaxFeF)) * 1.e3 ! = mmol/m3 
2679   !! 
2680   !! This assumes that the "excess" iron is dissipated with a time 
2681   !! scale of 1 day; seems reasonable to me ... (famous last words) 
2682   !! 
2683   ffescav = ffescav + fdeltaFe ! = mmol/m3/d 
2684   !! 
2685   # if defined key_deep_fe_fix 
2686   !! AXY (17/01/13) 
2687   !! stop scavenging for iron concentrations below 0.5 umol / m3 
2688   !! at depths greater than 1000 m; this aims to end MEDUSA's 
2689   !! continual loss of iron at depth without impacting things 
2690   !! at the surface too much; the justification for this is that 
2691   !! it appears to be what Mick Follows et al. do in their work 
2692   !! (as evidenced by the iron initial condition they supplied 
2693   !! me with); to be honest, it looks like Follow et al. do this 
2694   !! at shallower depths than 1000 m, but I'll stick with this 
2695   !! for now; I suspect that this seemingly arbitrary approach 
2696   !! effectively "parameterises" the particlebased scavenging 
2697   !! rates that other models use (i.e. at depth there are no 
2698   !! sinking particles, so scavenging stops); it might be fun 
2699   !! justifying this in a paper though! 
2700   !! 
2701   if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then 
2702   ffescav = 0. 
2703   endif 
2704   # endif 
2705   !! 
2706   elseif (jiron.eq.2) then 
2707   !! 
2708   !! Scheme 2: Moore et al. (2004) 
2709   !! This scheme includes a single scavenging term that accounts for 
2710   !! both suspended and sinking particles in the water column; this 
2711   !! term scavenges total iron rather than "free" iron 
2712   !! 
2713   !! 
2714   !! total iron concentration (mmol Fe / m3 > umol Fe / m3) 
2715   xFeT = zfer * 1.e3 
2716   !! 
2717   !! this has a base scavenging rate (12% / y) which is modified by local 
2718   !! particle concentration and sinking flux (and dust  but I'm ignoring 
2719   !! that here for now) and which is accelerated when Fe concentration gets 
2720   !! 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as concentrations 
2721   !! below 0.4 nM (= 0.4 umol/m3 = 0.0004 mmol/m3) 
2722   !! 
2723   !! base scavenging rate (0.12 / y) 
2724   fbase_scav = 0.12 / 365.25 
2725   !! 
2726   !! calculate sinking particle part of scaling factor 
2727   !! this takes local fast sinking carbon (mmol C / m2 / d) and 
2728   !! gets it into nmol C / cm3 / s ("rdt" below is the number of seconds in 
2729   !! a model timestep) 
2730   !! 
2731   !! fscal_sink = ffastc(ji,jj) * 1.e2 / (86400.) 
2732   !! 
2733   !! ... actually, rereading Moore et al.'s equations, it looks like he uses 
2734   !! his sinking flux directly, without scaling it by timestep or anything, 
2735   !! so I'll copy this here ... 
2736   !! 
2737   fscal_sink = ffastc(ji,jj) * 1.e2 
2738   !! 
2739   !! calculate particle part of scaling factor 
2740   !! this totals up the carbon in suspended particles (Pn, Pd, Zmi, Zme, D), 
2741   !! which comes out in mmol C / m3 (= nmol C / cm3), and then multiplies it 
2742   !! by a magic factor, 0.002, to get it into nmol C / cm2 / s 
2743   !! 
2744   fscal_part = ((xthetapn * zphn) + (xthetapd * zphd) + (xthetazmi * zzmi) + & 
2745   (xthetazme * zzme) + (xthetad * zdet)) * 0.002 
2746   !! 
2747   !! calculate scaling factor for base scavenging rate 
2748   !! this uses the (now correctly scaled) sinking flux and standing 
2749   !! particle concentration, divides through by some sort of reference 
2750   !! value (= 0.0066 nmol C / cm2 / s) and then uses this, or not if its 
2751   !! too high, to rescale the base scavenging rate 
2752   !! 
2753   fscal_scav = fbase_scav * min(((fscal_sink + fscal_part) / 0.0066), 4.0) 
2754   !! 
2755   !! the resulting scavenging rate is then scaled further according to the 
2756   !! local iron concentration (i.e. diminished in low iron regions; enhanced 
2757   !! in high iron regions; less alone in intermediate iron regions) 
2758   !! 
2759   if (xFeT.lt.0.4) then 
2760   !! 
2761   !! low iron region 
2762   !! 
2763   fscal_scav = fscal_scav * (xFeT / 0.4) 
2764   !! 
2765   elseif (xFeT.gt.0.6) then 
2766   !! 
2767   !! high iron region 
2768   !! 
2769   fscal_scav = fscal_scav + ((xFeT / 0.6) * (6.0 / 1.4)) 
2770   !! 
2771   else 
2772   !! 
2773   !! intermediate iron region: do nothing 
2774   !! 
2775   endif 
2776   !! 
2777   !! apply the calculated scavenging rate ... 
2778   !! 
2779   ffescav = fscal_scav * zfer 
2780   !! 
2781   elseif (jiron.eq.3) then 
2782   !! 
2783   !! Scheme 3: Moore et al. (2008) 
2784   !! This scheme includes a single scavenging term that accounts for 
2785   !! sinking particles in the water column, and includes organic C, 
2786   !! biogenic opal, calcium carbonate and dust in this (though the 
2787   !! latter is ignored here until I work out what units the incoming 
2788   !! "dust" flux is in); this term scavenges total iron rather than 
2789   !! "free" iron 
2790   !! 
2791   !! 
2792   !! total iron concentration (mmol Fe / m3 > umol Fe / m3) 
2793   xFeT = zfer * 1.e3 
2794   !! 
2795   !! this has a base scavenging rate which is modified by local 
2796   !! particle sinking flux (including dust  but I'm ignoring that 
2797   !! here for now) and which is accelerated when Fe concentration 
2798   !! is > 0.6 nM (= 0.6 umol/m3 = 0.0006 mmol/m3), and decreased as 
2799   !! concentrations < 0.5 nM (= 0.5 umol/m3 = 0.0005 mmol/m3) 
2800   !! 
2801   !! base scavenging rate (Fe_b in paper; units may be wrong there) 
2802   fbase_scav = 0.00384 ! (ng)^1 cm 
2803   !! 
2804   !! calculate sinking particle part of scaling factor; this converts 
2805   !! mmol / m2 / d fluxes of organic carbon, silicon and calcium 
2806   !! carbonate into ng / cm2 / s fluxes; it is assumed here that the 
2807   !! mass conversions simply consider the mass of the main element 
2808   !! (C, Si and Ca) and *not* the mass of the molecules that they are 
2809   !! part of; Moore et al. (2008) is unclear on the conversion that 
2810   !! should be used 
2811   !! 
2812   !! milli > nano; mol > gram; /m2 > /cm2; /d > /s 
2813   fscal_csink = (ffastc(ji,jj) * 1.e6 * xmassc * 1.e4 / 86400.) ! ng C / cm2 / s 
2814   fscal_sisink = (ffastsi(ji,jj) * 1.e6 * xmasssi * 1.e4 / 86400.) ! ng Si / cm2 / s 
2815   fscal_casink = (ffastca(ji,jj) * 1.e6 * xmassca * 1.e4 / 86400.) ! ng Ca / cm2 / s 
2816   !! 
2817   !! sum up these sinking fluxes and convert to ng / cm by dividing 
2818   !! through by a sinking rate of 100 m / d = 1.157 cm / s 
2819   fscal_sink = ((fscal_csink * 6.) + fscal_sisink + fscal_casink) / & 
2820   (100. * 1.e3 / 86400) ! ng / cm 
2821   !! 
2822   !! now calculate the scavenging rate based upon the base rate and 
2823   !! this particle flux scaling; according to the published units, 
2824   !! the result actually has *no* units, but as it must be expressed 
2825   !! per unit time for it to make any sense, I'm assuming a missing 
2826   !! "per second" 
2827   fscal_scav = fbase_scav * fscal_sink ! / s 
2828   !! 
2829   !! the resulting scavenging rate is then scaled further according to the 
2830   !! local iron concentration (i.e. diminished in low iron regions; enhanced 
2831   !! in high iron regions; less alone in intermediate iron regions) 
2832   !! 
2833   if (xFeT.lt.0.5) then 
2834   !! 
2835   !! low iron region (0.5 instead of the 0.4 in Moore et al., 2004) 
2836   !! 
2837   fscal_scav = fscal_scav * (xFeT / 0.5) 
2838   !! 
2839   elseif (xFeT.gt.0.6) then 
2840   !! 
2841   !! high iron region (functional form different in Moore et al., 2004) 
2842   !! 
2843   fscal_scav = fscal_scav + ((xFeT  0.6) * 0.00904) 
2844   !! 
2845   else 
2846   !! 
2847   !! intermediate iron region: do nothing 
2848   !! 
2849   endif 
2850   !! 
2851   !! apply the calculated scavenging rate ... 
2852   !! 
2853   ffescav = fscal_scav * zfer 
2854   !! 
2855   elseif (jiron.eq.4) then 
2856   !! 
2857   !! Scheme 4: Galbraith et al. (2010) 
2858   !! This scheme includes two scavenging terms, one for organic, 
2859   !! particlebased scavenging, and another for inorganic scavenging; 
2860   !! both terms scavenge "free" iron only 
2861   !! 
2862   !! 
2863   !! Galbraith et al. (2010) present a more straightforward outline of 
2864   !! the scheme in Parekh et al. (2005) ... 
2865   !! 
2866   !! sinking particulate carbon available for scavenging 
2867   !! this assumes a sinking rate of 100 m / d (Moore & Braucher, 2008), 
2868   xCscav1 = (ffastc(ji,jj) * xmassc) / 100. ! = mg C / m3 
2869   !! 
2870   !! scale by Honeyman et al. (1981) exponent coefficient 
2871   !! multiply by 1.e3 to express C flux in g C rather than mg C 
2872   xCscav2 = (xCscav1 * 1.e3)**0.58 
2873   !! 
2874   !! multiply by Galbraith et al. (2010) scavenging rate 
2875   xk_org = 0.5 ! ((g C m/3)^1) / d 
2876   xORGscav = xk_org * xCscav2 * xFeF 
2877   !! 
2878   !! Galbraith et al. (2010) also include an inorganic bit ... 
2879   !! 
2880   !! this occurs at a fixed rate, again based on the availability of 
2881   !! "free" iron 
2882   !! 
2883   !! k_inorg = 1000 d**1 nmol Fe**0.5 kg**0.5 
2884   !! 
2885   !! to implement this here, scale xFeF by 1026 to put in units of 
2886   !! umol/m3 which approximately equal nmol/kg 
2887   !! 
2888   xk_inorg = 1000. ! ((nmol Fe / kg)^1.5) 
2889   xINORGscav = (xk_inorg * (xFeF * 1026.)**1.5) * 1.e3 
2890   !! 
2891   !! sum these two terms together 
2892   ffescav = xORGscav + xINORGscav 
2893   else 
2894   !! 
2895   !! No Scheme: you coward! 
2896   !! This scheme puts its head in the sand and eskews any decision about 
2897   !! how iron is removed from the ocean; prepare to get deluged in iron 
2898   !! you fool! 
2899   !! 
2900   ffescav = 0. 
2901   endif 
2902   
2903   !! 
2904   !! Other iron cycle processes 
2905   !! 
2906   !! 
2907   !! aeolian iron deposition 
2908   if (jk.eq.1) then 
2909   !! zirondep is in mmolFe / m2 / day 
2910   !! ffetop is in mmoldissolvedFe / m3 / day 
2911   ffetop = zirondep(ji,jj) * xfe_sol / fthk 
2912   else 
2913   ffetop = 0.0 
2914   endif 
2915   !! 
2916   !! seafloor iron addition 
2917   !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down 
2918   !! if (jk.eq.(mbathy(ji,jj)1).AND.jk.lt.i1100) then 
2919   if ((jk.eq.jmbathy).AND.jk.le.i0500) then 
2920   !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a 
2921   !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value 
2922   !! but apply it everywhere 
2923   !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 137) 
2924   ffebot = (xfe_sed / fthk) 
2925   else 
2926   ffebot = 0.0 
2927   endif 
2928   
2929   !! AXY (16/12/09): remove iron addition/removal processes 
2930   !! For the purposes of the quarter degree run, the iron cycle is being pegged to the 
2931   !! initial condition supplied by Mick Follows via restoration with a 30 day period; 
2932   !! iron addition at the seafloor is still permitted with the idea that this extra 
2933   !! iron will be removed by the restoration away from the source 
2934   !! ffescav = 0.0 
2935   !! ffetop = 0.0 
2936   !! ffebot = 0.0 
2937   
2938   # if defined key_debug_medusa 
2939   !! report miscellaneous calculations 
2940   if (idf.eq.1.AND.idfval.eq.1) then 
2941   IF (lwp) write (numout,*) '' 
2942   IF (lwp) write (numout,*) 'xfe_sol = ', xfe_sol 
2943   IF (lwp) write (numout,*) 'xfe_mass = ', xfe_mass 
2944   IF (lwp) write (numout,*) 'ffetop(',jk,') = ', ffetop 
2945   IF (lwp) write (numout,*) 'ffebot(',jk,') = ', ffebot 
2946   IF (lwp) write (numout,*) 'xFree(',jk,') = ', xFree(ji,jj) 
2947   IF (lwp) write (numout,*) 'ffescav(',jk,') = ', ffescav 
2948   endif 
2949   # endif 
2950   
2951   !! 
2952   !! Miscellaneous 
2953   !! 
2954   !! 
2955   !! diatom frustule dissolution 
2956   fsdiss = xsdiss * zpds 
2957   
2958   # if defined key_debug_medusa 
2959   !! report miscellaneous calculations 
2960   if (idf.eq.1.AND.idfval.eq.1) then 
2961   IF (lwp) write (numout,*) '' 
2962   IF (lwp) write (numout,*) 'fsdiss(',jk,') = ', fsdiss 
2963   endif 
2964   # endif 
2965   
2966   !! 
2967   !! Slow detritus creation 
2968   !! 
2969   !! this variable integrates the creation of slow sinking detritus 
2970   !! to allow the split between fast and slow detritus to be 
2971   !! diagnosed 
2972   fslown = fdpn + fdzmi + ((1.0  xfdfrac1) * fdpd) + & 
2973   ((1.0  xfdfrac2) * fdzme) + ((1.0  xbetan) * (finmi + finme)) 
2974   !! 
2975   !! this variable records the slow detrital sinking flux at this 
2976   !! particular depth; it is used in the output of this flux at 
2977   !! standard depths in the diagnostic outputs; needs to be 
2978   !! adjusted from per second to per day because of parameter vsed 
2979   fslownflux(ji,jj) = zdet * vsed * 86400. 
2980   # if defined key_roam 
2981   !! 
2982   !! and the same for detrital carbon 
2983   fslowc = (xthetapn * fdpn) + (xthetazmi * fdzmi) + & 
2984   (xthetapd * (1.0  xfdfrac1) * fdpd) + & 
2985   (xthetazme * (1.0  xfdfrac2) * fdzme) + & 
2986   ((1.0  xbetac) * (ficmi + ficme)) 
2987   !! 
2988   !! this variable records the slow detrital sinking flux at this 
2989   !! particular depth; it is used in the output of this flux at 
2990   !! standard depths in the diagnostic outputs; needs to be 
2991   !! adjusted from per second to per day because of parameter vsed 
2992   fslowcflux(ji,jj) = zdtc * vsed * 86400. 
2993   # endif 
2994   
2995   !! 
2996   !! Nutrient regeneration 
2997   !! this variable integrates total nitrogen regeneration down the 
2998   !! watercolumn; its value is stored and output as a 2D diagnostic; 
2999   !! the corresponding dissolution flux of silicon (from sources 
3000   !! other than fast detritus) is also integrated; note that, 
3001   !! confusingly, the linear loss terms from plankton compartments 
3002   !! are labelled as fdX2 when one might have expected fdX or fdX1 
3003   !! 
3004   !! 
3005   !! nitrogen 
3006   fregen = (( (xphi * (fgmipn + fgmid)) + & ! messy feeding 
3007   (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + & ! messy feeding 
3008   fmiexcr + fmeexcr + fdd + & ! excretion + D remin. 
3009   fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk) ! linear mortality 
3010   !! 
3011   !! silicon 
3012   fregensi = (( fsdiss + ((1.0  xfdfrac1) * fdpds) + & ! dissolution + nonlin. mortality 
3013   ((1.0  xfdfrac3) * fgmepds) + & ! egestion by zooplankton 
3014   fdpds2) * fthk) ! linear mortality 
3015   # if defined key_roam 
3016   !! 
3017   !! carbon 
3018   fregenc = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) + & ! messy feeding 
3019   (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & ! messy feeding 
3020   (xthetazmi * fgmezmi) + fgmedc)) + & ! messy feeding 
3021   fmiresp + fmeresp + fddc + & ! respiration + D remin. 
3022   (xthetapn * fdpn2) + (xthetapd * fdpd2) + & ! linear mortality 
3023   (xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk) ! linear mortality 
3024   # endif 
3025   
3026   !! 
3027   !! Fastsinking detritus terms 
3028   !! "local" variables declared so that conservation can be checked; 
3029   !! the calculated terms are added to the fastsinking flux later on 
3030   !! only after the flux entering this level has experienced some 
3031   !! remineralisation 
3032   !! note: these fluxes need to be scaled by the level thickness 
3033   !! 
3034   !! 
3035   !! nitrogen: diatom and mesozooplankton mortality 
3036   ftempn = b0 * ((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) 
3037   !! 
3038   !! silicon: diatom mortality and grazed diatoms 
3039   ftempsi = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds)) 
3040   !! 
3041   !! iron: diatom and mesozooplankton mortality 
3042   ftempfe = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn) 
3043   !! 
3044   !! carbon: diatom and mesozooplankton mortality 
3045   ftempc = b0 * ((xfdfrac1 * xthetapd * fdpd) + & 
3046   (xfdfrac2 * xthetazme * fdzme)) 
3047   !! 
3048   # if defined key_roam 
3049   if (jrratio.eq.0) then 
3050   !! CaCO3: latitudinallybased fraction of total primary production 
3051   !! absolute latitude of current grid cell 
3052   flat = abs(gphit(ji,jj)) 
3053   !! 0.10 at equator; 0.02 at pole 
3054   fcaco3 = xcaco3a + ((xcaco3b  xcaco3a) * ((90.0  flat) / 90.0)) 
3055   elseif (jrratio.eq.1) then 
3056   !! CaCO3: Ridgwell et al. (2007) submodel, version 1 
3057   !! this uses SURFACE omega calcite to regulate rain ratio 
3058   if (f_omcal(ji,jj).ge.1.0) then 
3059   fq1 = (f_omcal(ji,jj)  1.0)**0.81 
3060   else 
3061   fq1 = 0. 
3062   endif 
3063   fcaco3 = xridg_r0 * fq1 
3064   elseif (jrratio.eq.2) then 
3065   !! CaCO3: Ridgwell et al. (2007) submodel, version 2 
3066   !! this uses FULL 3D omega calcite to regulate rain ratio 
3067   if (f3_omcal(ji,jj,jk).ge.1.0) then 
3068   fq1 = (f3_omcal(ji,jj,jk)  1.0)**0.81 
3069   else 
3070   fq1 = 0. 
3071   endif 
3072   fcaco3 = xridg_r0 * fq1 
3073   endif 
3074   # else 
3075   !! CaCO3: latitudinallybased fraction of total primary production 
3076   !! absolute latitude of current grid cell 
3077   flat = abs(gphit(ji,jj)) 
3078   !! 0.10 at equator; 0.02 at pole 
3079   fcaco3 = xcaco3a + ((xcaco3b  xcaco3a) * ((90.0  flat) / 90.0)) 
3080   # endif 
3081   !! AXY (09/03/09): convert CaCO3 production from function of 
3082   !! primary production into a function of fastsinking material; 
3083   !! technically, this is what Dunne et al. (2007) do anyway; they 
3084   !! convert total primary production estimated from surface 
3085   !! chlorophyll to an export flux for which they apply conversion 
3086   !! factors to estimate the various elemental fractions (Si, Ca) 
3087   ftempca = ftempc * fcaco3 
3088   
3089   # if defined key_debug_medusa 
3090   !! integrate total fast detritus production 
3091   if (idf.eq.1) then 
3092   fifd_n(ji,jj) = fifd_n(ji,jj) + (ftempn * fthk) 
3093   fifd_si(ji,jj) = fifd_si(ji,jj) + (ftempsi * fthk) 
3094   fifd_fe(ji,jj) = fifd_fe(ji,jj) + (ftempfe * fthk) 
3095   # if defined key_roam 
3096   fifd_c(ji,jj) = fifd_c(ji,jj) + (ftempc * fthk) 
3097   # endif 
3098   endif 
3099   
3100   !! report quantities of fastsinking detritus for each component 
3101   if (idf.eq.1.AND.idfval.eq.1) then 
3102   IF (lwp) write (numout,*) '' 
3103   IF (lwp) write (numout,*) 'fdpd(',jk,') = ', fdpd 
3104   IF (lwp) write (numout,*) 'fdzme(',jk,') = ', fdzme 
3105   IF (lwp) write (numout,*) 'ftempn(',jk,') = ', ftempn 
3106   IF (lwp) write (numout,*) 'ftempsi(',jk,') = ', ftempsi 
3107   IF (lwp) write (numout,*) 'ftempfe(',jk,') = ', ftempfe 
3108   IF (lwp) write (numout,*) 'ftempc(',jk,') = ', ftempc 
3109   IF (lwp) write (numout,*) 'ftempca(',jk,') = ', ftempca 
3110   IF (lwp) write (numout,*) 'flat(',jk,') = ', flat 
3111   IF (lwp) write (numout,*) 'fcaco3(',jk,') = ', fcaco3 
3112   endif 
3113   # endif 
3114   
3115   !! 
3116   !! This version of MEDUSA offers a choice of three methods for 
3117   !! handling the remineralisation of fast detritus. All three 
3118   !! do so in broadly the same way: 
3119   !! 
3120   !! 1. Fast detritus is stored as a 2D array [ ffastX ] 
3121   !! 2. Fast detritus is added levelbylevel [ ftempX ] 
3122   !! 3. Fast detritus is not remineralised in the top box [ freminX ] 
3123   !! 4. Remaining fast detritus is remineralised in the bottom [ fsedX ] 
3124   !! box 
3125   !! 
3126   !! The three remineralisation methods are: 
3127   !! 
3128   !! 1. Ballast model (i.e. that published in Yool et al., 2011) 
3129   !! (1b. Ballastsansballast model) 
3130   !! 2. Martin et al. (1987) 
3131   !! 3. Henson et al. (2011) 
3132   !! 
3133   !! The first of these couples C, N and Fe remineralisation to 
3134   !! the remineralisation of particulate Si and CaCO3, but the 
3135   !! latter two treat remineralisation of C, N, Fe, Si and CaCO3 
3136   !! completely separately. At present a switch within the code 
3137   !! regulates which submodel is used, but this should be moved 
3138   !! to the namelist file. 
3139   !! 
3140   !! The ballastsansballast submodel is an original development 
3141   !! feature of MEDUSA in which the ballast submodel's general 
3142   !! framework and parameterisation is used, but in which there 
3143   !! is no protection of organic material afforded by ballasting 
3144   !! minerals. While similar, it is not the same as the Martin 
3145   !! et al. (1987) submodel. 
3146   !! 
3147   !! Since the three submodels behave the same in terms of 
3148   !! accumulating sinking material and remineralising it all at 
3149   !! the seafloor, these portions of the code below are common to 
3150   !! all three. 
3151   !! 
3152   
3153   if (jexport.eq.1) then 
3154   !!====================================================================== 
3155   !! BALLAST SUBMODEL 
3156   !!====================================================================== 
3157   !! 
3158   !! 
3159   !! Fastsinking detritus fluxes, pt. 1: REMINERALISATION 
3160   !! aside from explicitly modelled, slowsinking detritus, the 
3161   !! model includes an implicit representation of detrital 
3162   !! particles that sink too quickly to be modelled with 
3163   !! explicit state variables; this sinking flux is instead 
3164   !! instantaneously remineralised down the water column using 
3165   !! the version of Armstrong et al. (2002)'s ballast model 
3166   !! used by Dunne et al. (2007); the version of this model 
3167   !! here considers silicon and calcium carbonate ballast 
3168   !! minerals; this section of the code redistributes the fast 
3169   !! sinking material generated locally down the water column; 
3170   !! this differs from Dunne et al. (2007) in that fast sinking 
3171   !! material is distributed at *every* level below that it is 
3172   !! generated, rather than at every level below some fixed 
3173   !! depth; this scheme is also different in that sinking material 
3174   !! generated in one level is aggregated with that generated by 
3175   !! shallower levels; this should make the ballast model more 
3176   !! selfconsistent (famous last words) 
3177   !! 
3178   !! 
3179   if (jk.eq.1) then 
3180   !! this is the SURFACE OCEAN BOX (no remineralisation) 
3181   !! 
3182   freminc = 0.0 
3183   freminn = 0.0 
3184   freminfe = 0.0 
3185   freminsi = 0.0 
3186   freminca = 0.0 
3187   elseif (jk.le.jmbathy) then 
3188   !! this is an OCEAN BOX (remineralise some material) 
3189   !! 
3190   !! set up CCD depth to be used depending on user choice 
3191   if (jocalccd.eq.0) then 
3192   !! use default CCD field 
3193   fccd_dep = ocal_ccd(ji,jj) 
3194   elseif (jocalccd.eq.1) then 
3195   !! use calculated CCD field 
3196   fccd_dep = f2_ccd_cal(ji,jj) 
3197   endif 
3198   !! 
3199   !! === organic carbon === 
3200   fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) 
3201   if (iball.eq.1) then 
3202   fq1 = (fq0 * xmassc) !! how much it weighs (mass) 
3203   fq2 = (ffastca(ji,jj) * xmassca) !! how much CaCO3 enters this box (mass) 
3204   fq3 = (ffastsi(ji,jj) * xmasssi) !! how much opal enters this box (mass) 
3205   fq4 = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C (mass) 
3206   !! this next term is calculated for C but used for N and Fe as well 
3207   !! it needs to be protected in case ALL C is protected 
3208   if (fq4.lt.fq1) then 
3209   fprotf = (fq4 / (fq1 + tiny(fq1))) !! protected fraction of total organic C (nondim) 
3210   else 
3211   fprotf = 1.0 !! all organic C is protected (nondim) 
3212   endif 
3213   fq5 = (1.0  fprotf) !! unprotected fraction of total organic C (nondim) 
3214   fq6 = (fq0 * fq5) !! how much organic C is unprotected (mol) 
3215   fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected C leaves this box (mol) 
3216   fq8 = (fq7 + (fq0 * fprotf)) !! how much total C leaves this box (mol) 
3217   freminc = (fq0  fq8) / fthk !! C remineralisation in this box (mol) 
3218   ffastc(ji,jj) = fq8 
3219   # if defined key_debug_medusa 
3220   !! report in/out/remin fluxes of carbon for this level 
3221   if (idf.eq.1.AND.idfval.eq.1) then 
3222   IF (lwp) write (numout,*) '' 
3223   IF (lwp) write (numout,*) 'totalC(',jk,') = ', fq1 
3224   IF (lwp) write (numout,*) 'prtctC(',jk,') = ', fq4 
3225   IF (lwp) write (numout,*) 'fprotf(',jk,') = ', fprotf 
3226   IF (lwp) write (numout,*) '' 
3227   IF (lwp) write (numout,*) 'IN C(',jk,') = ', fq0 
3228   IF (lwp) write (numout,*) 'LOST C(',jk,') = ', freminc * fthk 
3229   IF (lwp) write (numout,*) 'OUT C(',jk,') = ', fq8 
3230   IF (lwp) write (numout,*) 'NEW C(',jk,') = ', ftempc * fthk 
3231   endif 
3232   # endif 
3233   else 
3234   fq1 = fq0 * exp((fthk / xfastc)) !! how much organic C leaves this box (mol) 
3235   freminc = (fq0  fq1) / fthk !! C remineralisation in this box (mol) 
3236   ffastc(ji,jj) = fq1 
3237   endif 
3238   !! 
3239   !! === organic nitrogen === 
3240   fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) 
3241   if (iball.eq.1) then 
3242   fq5 = (1.0  fprotf) !! unprotected fraction of total organic N (nondim) 
3243   fq6 = (fq0 * fq5) !! how much organic N is unprotected (mol) 
3244   fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected N leaves this box (mol) 
3245   fq8 = (fq7 + (fq0 * fprotf)) !! how much total N leaves this box (mol) 
3246   freminn = (fq0  fq8) / fthk !! N remineralisation in this box (mol) 
3247   ffastn(ji,jj) = fq8 
3248   # if defined key_debug_medusa 
3249   !! report in/out/remin fluxes of carbon for this level 
3250   if (idf.eq.1.AND.idfval.eq.1) then 
3251   IF (lwp) write (numout,*) '' 
3252   IF (lwp) write (numout,*) 'totalN(',jk,') = ', fq1 
3253   IF (lwp) write (numout,*) 'prtctN(',jk,') = ', fq4 
3254   IF (lwp) write (numout,*) 'fprotf(',jk,') = ', fprotf 
3255   IF (lwp) write (numout,*) '' 
3256   if (freminn < 0.0) then 
3257   IF (lwp) write (numout,*) '** FREMIN ERROR **' 
3258   endif 
3259   IF (lwp) write (numout,*) 'IN N(',jk,') = ', fq0 
3260   IF (lwp) write (numout,*) 'LOST N(',jk,') = ', freminn * fthk 
3261   IF (lwp) write (numout,*) 'OUT N(',jk,') = ', fq8 
3262   IF (lwp) write (numout,*) 'NEW N(',jk,') = ', ftempn * fthk 
3263   endif 
3264   # endif 
3265   else 
3266   fq1 = fq0 * exp((fthk / xfastc)) !! how much organic N leaves this box (mol) 
3267   freminn = (fq0  fq1) / fthk !! N remineralisation in this box (mol) 
3268   ffastn(ji,jj) = fq1 
3269   endif 
3270   !! 
3271   !! === organic iron === 
3272   fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) 
3273   if (iball.eq.1) then 
3274   fq5 = (1.0  fprotf) !! unprotected fraction of total organic Fe (nondim) 
3275   fq6 = (fq0 * fq5) !! how much organic Fe is unprotected (mol) 
3276   fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected Fe leaves this box (mol) 
3277   fq8 = (fq7 + (fq0 * fprotf)) !! how much total Fe leaves this box (mol) 
3278   freminfe = (fq0  fq8) / fthk !! Fe remineralisation in this box (mol) 
3279   ffastfe(ji,jj) = fq8 
3280   else 
3281   fq1 = fq0 * exp((fthk / xfastc)) !! how much total Fe leaves this box (mol) 
3282   freminfe = (fq0  fq1) / fthk !! Fe remineralisation in this box (mol) 
3283   ffastfe(ji,jj) = fq1 
3284   endif 
3285   !! 
3286   !! === biogenic silicon === 
3287   fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) 
3288   fq1 = fq0 * exp((fthk / xfastsi)) !! how much opal leaves this box (mol) 
3289   freminsi = (fq0  fq1) / fthk !! Si remineralisation in this box (mol) 
3290   ffastsi(ji,jj) = fq1 
3291   !! 
3292   !! === biogenic calcium carbonate === 
3293   fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) 
3294   if (fdep.le.fccd_dep) then 
3295   !! whole grid cell above CCD 
3296   fq1 = fq0 !! above lysocline, no Ca dissolves (mol) 
3297   freminca = 0.0 !! above lysocline, no Ca dissolves (mol) 
3298   fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) 
3299   elseif (fdep.ge.fccd_dep) then 
3300   !! whole grid cell below CCD 
3301   fq1 = fq0 * exp((fthk / xfastca)) !! how much CaCO3 leaves this box (mol) 
3302   freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
3303   else 
3304   !! partial grid cell below CCD 
3305   fq2 = fdep1  fccd_dep !! amount of grid cell below CCD (m) 
3306   fq1 = fq0 * exp((fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) 
3307   freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
3308   endif 
3309   ffastca(ji,jj) = fq1 
3310   else 
3311   !! this is BELOW THE LAST OCEAN BOX (do nothing) 
3312   freminc = 0.0 
3313   freminn = 0.0 
3314   freminfe = 0.0 
3315   freminsi = 0.0 
3316   freminca = 0.0 
3317   endif 
3318   
3319   elseif (jexport.eq.2.or.jexport.eq.3) then 
3320   if (jexport.eq.2) then 
3321   !!====================================================================== 
3322   !! MARTIN ET AL. (1987) SUBMODEL 
3323   !!====================================================================== 
3324   !! 
3325   !! 
3326   !! This submodel uses the classic Martin et al. (1987) curve 
3327   !! to determine the attenuation of fastsinking detritus down 
3328   !! the water column. All three organic elements, C, N and Fe, 
3329   !! are handled identically, and their quantities in sinking 
3330   !! particles attenuate according to a power relationship 
3331   !! governed by parameter "b". This is assigned a canonical 
3332   !! value of 0.858. Biogenic opal and calcium carbonate are 
3333   !! attentuated using the same function as in the ballast 
3334   !! submodel 
3335   !! 
3336   !! 
3337   fb_val = 0.858 
3338   elseif (jexport.eq.3) then 
3339   !!====================================================================== 
3340   !! HENSON ET AL. (2011) SUBMODEL 
3341   !!====================================================================== 
3342   !! 
3343   !! 
3344   !! This submodel reconfigures the Martin et al. (1987) curve by 
3345   !! allowing the "b" value to vary geographically. Its value is 
3346   !! set, following Henson et al. (2011), as a function of local 
3347   !! sea surface temperature: 
3348   !! b = 1.06 + (0.024 * SST) 
3349   !! This means that remineralisation length scales are longer in 
3350   !! warm, tropical areas and shorter in cold, polar areas. This 
3351   !! does seem backtofront (i.e. one would expect GREATER 
3352   !! remineralisation in warmer waters), but is an outcome of 
3353   !! analysis of sediment trap data, and it may reflect details 
3354   !! of ecosystem structure that pertain to particle production 
3355   !! rather than simply Q10. 
3356   !! 
3357   !! 
3358   fl_sst = tsn(ji,jj,1,jp_tem) 
3359   fb_val = 1.06 + (0.024 * fl_sst) 
3360   endif 
3361   !! 
3362   if (jk.eq.1) then 
3363   !! this is the SURFACE OCEAN BOX (no remineralisation) 
3364   !! 
3365   freminc = 0.0 
3366   freminn = 0.0 
3367   freminfe = 0.0 
3368   freminsi = 0.0 
3369   freminca = 0.0 
3370   elseif (jk.le.jmbathy) then 
3371   !! this is an OCEAN BOX (remineralise some material) 
3372   !! 
3373   !! === organic carbon === 
3374   fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) 
3375   fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic C leaves this box (mol) 
3376   freminc = (fq0  fq1) / fthk !! C remineralisation in this box (mol) 
3377   ffastc(ji,jj) = fq1 
3378   !! 
3379   !! === organic nitrogen === 
3380   fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) 
3381   fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic N leaves this box (mol) 
3382   freminn = (fq0  fq1) / fthk !! N remineralisation in this box (mol) 
3383   ffastn(ji,jj) = fq1 
3384   !! 
3385   !! === organic iron === 
3386   fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) 
3387   fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic Fe leaves this box (mol) 
3388   freminfe = (fq0  fq1) / fthk !! Fe remineralisation in this box (mol) 
3389   ffastfe(ji,jj) = fq1 
3390   !! 
3391   !! === biogenic silicon === 
3392   fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) 
3393   fq1 = fq0 * exp((fthk / xfastsi)) !! how much opal leaves this box (mol) 
3394   freminsi = (fq0  fq1) / fthk !! Si remineralisation in this box (mol) 
3395   ffastsi(ji,jj) = fq1 
3396   !! 
3397   !! === biogenic calcium carbonate === 
3398   fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) 
3399   if (fdep.le.ocal_ccd(ji,jj)) then 
3400   !! whole grid cell above CCD 
3401   fq1 = fq0 !! above lysocline, no Ca dissolves (mol) 
3402   freminca = 0.0 !! above lysocline, no Ca dissolves (mol) 
3403   fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) 
3404   elseif (fdep.ge.ocal_ccd(ji,jj)) then 
3405   !! whole grid cell below CCD 
3406   fq1 = fq0 * exp((fthk / xfastca)) !! how much CaCO3 leaves this box (mol) 
3407   freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
3408   else 
3409   !! partial grid cell below CCD 
3410   fq2 = fdep1  ocal_ccd(ji,jj) !! amount of grid cell below CCD (m) 
3411   fq1 = fq0 * exp((fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) 
3412   freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
3413   endif 
3414   ffastca(ji,jj) = fq1 
3415   else 
3416   !! this is BELOW THE LAST OCEAN BOX (do nothing) 
3417   freminc = 0.0 
3418   freminn = 0.0 
3419   freminfe = 0.0 
3420   freminsi = 0.0 
3421   freminca = 0.0 
3422   endif 
3423   
3424   endif 
3425   
3426   !! 
3427   !! Fastsinking detritus fluxes, pt. 2: UPDATE FAST FLUXES 
3428   !! here locally calculated additions to the fastsinking flux are added 
3429   !! to the total fastsinking flux; this is done here such that material 
3430   !! produced in a particular layer is only remineralised below this 
3431   !! layer 
3432   !! 
3433   !! 
3434   !! add sinking material generated in this layer to running totals 
3435   !! 
3436   !! === organic carbon === (diatom and mesozooplankton mortality) 
3437   ffastc(ji,jj) = ffastc(ji,jj) + (ftempc * fthk) 
3438   !! 
3439   !! === organic nitrogen === (diatom and mesozooplankton mortality) 
3440   ffastn(ji,jj) = ffastn(ji,jj) + (ftempn * fthk) 
3441   !! 
3442   !! === organic iron === (diatom and mesozooplankton mortality) 
3443   ffastfe(ji,jj) = ffastfe(ji,jj) + (ftempfe * fthk) 
3444   !! 
3445   !! === biogenic silicon === (diatom mortality and grazed diatoms) 
3446   ffastsi(ji,jj) = ffastsi(ji,jj) + (ftempsi * fthk) 
3447   !! 
3448   !! === biogenic calcium carbonate === (latitudinallybased fraction of total primary production) 
3449   ffastca(ji,jj) = ffastca(ji,jj) + (ftempca * fthk) 
3450   
3451   !! 
3452   !! Fastsinking detritus fluxes, pt. 3: SEAFLOOR 
3453   !! remineralise all remaining fastsinking detritus to dissolved 
3454   !! nutrients; the sedimentation fluxes calculated here allow the 
3455   !! separation of what's remineralised sinking through the final 
3456   !! ocean box from that which is added to the final box by the 
3457   !! remineralisation of material that reaches the seafloor (i.e. 
3458   !! the model assumes that *all* material that hits the seafloor 
3459   !! is remineralised and that none is permanently buried; hey, 
3460   !! this is a giant GCM model that can't be run for long enough 
3461   !! to deal with burial fluxes!) 
3462   !! 
3463   !! in a change to this process, in part so that MEDUSA behaves 
3464   !! a little more like ERSEM et al., fastsinking detritus (N, Fe 
3465   !! and C) is converted to slow sinking detritus at the seafloor 
3466   !! instead of being remineralised; the rationale is that in 
3467   !! shallower shelf regions (... that are not fully mixed!) this 
3468   !! allows the detrital material to return slowly to dissolved 
3469   !! nutrient rather than instantaneously as now; the alternative 
3470   !! would be to explicitly handle seafloor organic material  a 
3471   !! headache I don't wish to experience at this point; note that 
3472   !! fastsinking Si and Ca detritus is just remineralised as 
3473   !! per usual 
3474   !! 
3475   !! AXY (13/01/12) 
3476   !! in a further change to this process, again so that MEDUSA is 
3477   !! a little more like ERSEM et al., material that reaches the 
3478   !! seafloor can now be added to sediment pools and stored for 
3479   !! slow release; there are new 2D arrays for organic nitrogen, 
3480   !! iron and carbon and inorganic silicon and carbon that allow 
3481   !! fast and slow detritus that reaches the seafloor to be held 
3482   !! and released back to the water column more slowly; these arrays 
3483   !! are transferred via the tracer restart files between repeat 
3484   !! submissions of the model 
3485   !! 
3486   !! 
3487   ffast2slowc = 0.0 
3488   ffast2slown = 0.0 
3489   ffast2slowfe = 0.0 
3490   !! 
3491   if (jk.eq.jmbathy) then 
3492   !! this is the BOTTOM OCEAN BOX (remineralise everything) 
3493   !! 
3494   !! AXY (17/01/12): tweaked to include benthos pools 
3495   !! 
3496   !! === organic carbon === 
3497   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
3498   freminc = freminc + (ffastc(ji,jj) / fthk) !! C remineralisation in this box (mol/m3) 
3499   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
3500   ffast2slowc = ffastc(ji,jj) / fthk !! fast C > slow C (mol/m3) 
3501   fslowc = fslowc + ffast2slowc 
3502   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
3503   f_fbenin_c(ji,jj) = ffastc(ji,jj) !! fast C > benthic C (mol/m2) 
3504   endif 
3505   fsedc(ji,jj) = ffastc(ji,jj) !! record seafloor C (mol/m2) 
3506   ffastc(ji,jj) = 0.0 
3507   !! 
3508   !! === organic nitrogen === 
3509   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
3510   freminn = freminn + (ffastn(ji,jj) / fthk) !! N remineralisation in this box (mol/m3) 
3511   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
3512   ffast2slown = ffastn(ji,jj) / fthk !! fast N > slow N (mol/m3) 
3513   fslown = fslown + ffast2slown 
3514   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
3515   f_fbenin_n(ji,jj) = ffastn(ji,jj) !! fast N > benthic N (mol/m2) 
3516   endif 
3517   fsedn(ji,jj) = ffastn(ji,jj) !! record seafloor N (mol/m2) 
3518   ffastn(ji,jj) = 0.0 
3519   !! 
3520   !! === organic iron === 
3521   if (jfdfate.eq.0 .and. jorgben.eq.0) then 
3522   freminfe = freminfe + (ffastfe(ji,jj) / fthk) !! Fe remineralisation in this box (mol/m3) 
3523   elseif (jfdfate.eq.1 .and. jorgben.eq.0) then 
3524   ffast2slowfe = ffastn(ji,jj) / fthk !! fast Fe > slow Fe (mol/m3) 
3525   elseif (jfdfate.eq.0 .and. jorgben.eq.1) then 
3526   f_fbenin_fe(ji,jj) = ffastfe(ji,jj) !! fast Fe > benthic Fe (mol/m2) 
3527   endif 
3528   fsedfe(ji,jj) = ffastfe(ji,jj) !! record seafloor Fe (mol/m2) 
3529   ffastfe(ji,jj) = 0.0 
3530   !! 
3531   !! === biogenic silicon === 
3532   if (jinorgben.eq.0) then 
3533   freminsi = freminsi + (ffastsi(ji,jj) / fthk) !! Si remineralisation in this box (mol/m3) 
3534   elseif (jinorgben.eq.1) then 
3535   f_fbenin_si(ji,jj) = ffastsi(ji,jj) !! fast Si > benthic Si (mol/m2) 
3536   endif 
3537   fsedsi(ji,jj) = ffastsi(ji,jj) !! record seafloor Si (mol/m2) 
3538   ffastsi(ji,jj) = 0.0 
3539   !! 
3540   !! === biogenic calcium carbonate === 
3541   if (jinorgben.eq.0) then 
3542   freminca = freminca + (ffastca(ji,jj) / fthk) !! Ca remineralisation in this box (mol/m3) 
3543   elseif (jinorgben.eq.1) then 
3544   f_fbenin_ca(ji,jj) = ffastca(ji,jj) !! fast Ca > benthic Ca (mol/m2) 
3545   endif 
3546   fsedca(ji,jj) = ffastca(ji,jj) !! record seafloor Ca (mol/m2) 
3547   ffastca(ji,jj) = 0.0 
3548   endif 
3549   
3550   # if defined key_debug_medusa 
3551   if (idf.eq.1) then 
3552   !! 
3553   !! Integrate total fast detritus remineralisation 
3554   !! 
3555   !! 
3556   fofd_n(ji,jj) = fofd_n(ji,jj) + (freminn * fthk) 
3557   fofd_si(ji,jj) = fofd_si(ji,jj) + (freminsi * fthk) 
3558   fofd_fe(ji,jj) = fofd_fe(ji,jj) + (freminfe * fthk) 
3559   # if defined key_roam 
3560   fofd_c(ji,jj) = fofd_c(ji,jj) + (freminc * fthk) 
3561   # endif 
3562   endif 
3563   # endif 
3564   
3565   !! 
3566   !! Sort out remineralisation tally of fastsinking detritus 
3567   !! 
3568   !! 
3569   !! update fastsinking regeneration arrays 
3570   fregenfast(ji,jj) = fregenfast(ji,jj) + (freminn * fthk) 
3571   fregenfastsi(ji,jj) = fregenfastsi(ji,jj) + (freminsi * fthk) 
3572   # if defined key_roam 
3573   fregenfastc(ji,jj) = fregenfastc(ji,jj) + (freminc * fthk) 
3574   # endif 
3575   
3576   !! 
3577   !! Benthic remineralisation fluxes 
3578   !! 
3579   !! 
3580   if (jk.eq.jmbathy) then 
3581   !! 
3582   !! organic components 
3583   if (jorgben.eq.1) then 
3584   f_benout_n(ji,jj) = xsedn * zn_sed_n(ji,jj) 
3585   f_benout_fe(ji,jj) = xsedfe * zn_sed_fe(ji,jj) 
3586   f_benout_c(ji,jj) = xsedc * zn_sed_c(ji,jj) 
3587   endif 
3588   !! 
3589   !! inorganic components 
3590   if (jinorgben.eq.1) then 
3591   f_benout_si(ji,jj) = xsedsi * zn_sed_si(ji,jj) 
3592   f_benout_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) 
3593   !! 
3594   !! account for CaCO3 that dissolves when it shouldn't 
3595   if ( fdep .le. fccd_dep ) then 
3596   f_benout_lyso_ca(ji,jj) = xsedca * zn_sed_ca(ji,jj) 
3597   endif 
3598   endif 
3599   endif 
3600   CALL flush(numout) 
3601   
3602   !!====================================================================== 
3603   !! LOCAL GRID CELL TRENDS 
3604   !!====================================================================== 
3605   !! 
3606   !! 
3607   !! Determination of trends 
3608   !! 
3609   !! 
3610   !! 
3611   !! chlorophyll 
3612   btra(jpchn) = b0 * ( & 
3613   + ((frn * fprn * zphn)  fgmipn  fgmepn  fdpn  fdpn2) * (fthetan / xxi) ) 
3614   btra(jpchd) = b0 * ( & 
3615   + ((frd * fprd * zphd)  fgmepd  fdpd  fdpd2) * (fthetad / xxi) ) 
3616   !! 
3617   !! 
3618   !! phytoplankton 
3619   btra(jpphn) = b0 * ( & 
3620   + (fprn * zphn)  fgmipn  fgmepn  fdpn  fdpn2 ) 
3621   btra(jpphd) = b0 * ( & 
3622   + (fprd * zphd)  fgmepd  fdpd  fdpd2 ) 
3623   btra(jppds) = b0 * ( & 
3624   + (fprds * zpds)  fgmepds  fdpds  fsdiss  fdpds2 ) 
3625   !! 
3626   !! 
3627   !! zooplankton 
3628   btra(jpzmi) = b0 * ( & 
3629   + fmigrow  fgmezmi  fdzmi  fdzmi2 ) 
3630   btra(jpzme) = b0 * ( & 
3631   + fmegrow  fdzme  fdzme2 ) 
3632   !! 
3633   !! 
3634   !! detritus 
3635   btra(jpdet) = b0 * ( & 
3636   + fdpn + ((1.0  xfdfrac1) * fdpd) & ! mort. losses 
3637   + fdzmi + ((1.0  xfdfrac2) * fdzme) & ! mort. losses 
3638   + ((1.0  xbetan) * (finmi + finme)) & ! assim. inefficiency 
3639    fgmid  fgmed  fdd & ! grazing and remin. 
3640   + ffast2slown ) ! seafloor fast>slow 
3641   !! 
3642   !! 
3643   !! dissolved inorganic nitrogen nutrient 
3644   fn_cons = 0.0 & 
3645    (fprn * zphn)  (fprd * zphd) ! primary production 
3646   fn_prod = 0.0 & 
3647   + (xphi * (fgmipn + fgmid)) & ! messy feeding remin. 
3648   + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) & ! messy feeding remin. 
3649   + fmiexcr + fmeexcr + fdd + freminn & ! excretion and remin. 
3650   + fdpn2 + fdpd2 + fdzmi2 + fdzme2 ! metab. losses 
3651   !! 
3652   !! riverine flux 
3653   if ( jriver_n .gt. 0 ) then 
3654   f_riv_loc_n = f_riv_n(ji,jj) * friver_dep(jk,jmbathy) / fthk 
3655   fn_prod = fn_prod + f_riv_loc_n 
3656   endif 
3657   !! 
3658   !! benthic remineralisation 
3659   if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then 
3660   fn_prod = fn_prod + (f_benout_n(ji,jj) / fthk) 
3661   endif 
3662   !! 
3663   btra(jpdin) = b0 * ( & 
3664   fn_prod + fn_cons ) 
3665   !! 
3666   fnit_cons(ji,jj) = fnit_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved nitrogen 
3667   fn_cons ) ) 
3668   fnit_prod(ji,jj) = fnit_prod(ji,jj) + ( fthk * ( & ! production of dissolved nitrogen 
3669   fn_prod ) ) 
3670   !! 
3671   !! 
3672   !! dissolved silicic acid nutrient 
3673   fs_cons = 0.0 & 
3674    (fprds * zpds) ! opal production 
3675   fs_prod = 0.0 & 
3676   + fsdiss & ! opal dissolution 
3677   + ((1.0  xfdfrac1) * fdpds) & ! mort. loss 
3678   + ((1.0  xfdfrac3) * fgmepds) & ! egestion of grazed Si 
3679   + freminsi + fdpds2 ! fast diss. and metab. losses 
3680   !! 
3681   !! riverine flux 
3682   if ( jriver_si .gt. 0 ) then 
3683   f_riv_loc_si = f_riv_si(ji,jj) * friver_dep(jk,jmbathy) / fthk 
3684   fs_prod = fs_prod + f_riv_loc_si 
3685   endif 
3686   !! 
3687   !! benthic remineralisation 
3688   if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then 
3689   fs_prod = fs_prod + (f_benout_si(ji,jj) / fthk) 
3690   endif 
3691   !! 
3692   btra(jpsil) = b0 * ( & 
3693   fs_prod + fs_cons ) 
3694   !! 
3695   fsil_cons(ji,jj) = fsil_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved silicon 
3696   fs_cons ) ) 
3697   fsil_prod(ji,jj) = fsil_prod(ji,jj) + ( fthk * ( & ! production of dissolved silicon 
3698   fs_prod ) ) 
3699   !! 
3700   !! 
3701   !! dissolved "iron" nutrient 
3702   btra(jpfer) = b0 * ( & 
3703   + (xrfn * btra(jpdin)) + ffetop + ffebot  ffescav ) 
3704   
3705   # if defined key_roam 
3706   !! 
3707   !! 
3708   !! AXY (26/11/08): implicit detrital carbon change 
3709   btra(jpdtc) = b0 * ( & 
3710   + (xthetapn * fdpn) + ((1.0  xfdfrac1) * (xthetapd * fdpd)) & ! mort. losses 
3711   + (xthetazmi * fdzmi) + ((1.0  xfdfrac2) * (xthetazme * fdzme)) & ! mort. losses 
3712   + ((1.0  xbetac) * (ficmi + ficme)) & ! assim. inefficiency 
3713    fgmidc  fgmedc  fddc & ! grazing and remin. 
3714   + ffast2slowc ) ! seafloor fast>slow 
3715   !! 
3716   !! 
3717   !! dissolved inorganic carbon 
3718   fc_cons = 0.0 & 
3719    (xthetapn * fprn * zphn)  (xthetapd * fprd * zphd) ! primary production 
3720   fc_prod = 0.0 & 
3721   + (xthetapn * xphi * fgmipn) + (xphi * fgmidc) & ! messy feeding remin 
3722   + (xthetapn * xphi * fgmepn) + (xthetapd * xphi * fgmepd) & ! messy feeding remin 
3723   + (xthetazmi * xphi * fgmezmi) + (xphi * fgmedc) & ! messy feeding remin 
3724   + fmiresp + fmeresp + fddc + freminc + (xthetapn * fdpn2) & ! resp., remin., losses 
3725   + (xthetapd * fdpd2) + (xthetazmi * fdzmi2) & ! losses 
3726   + (xthetazme * fdzme2) ! losses 
3727   !! 
3728   !! riverine flux 
3729   if ( jriver_c .gt. 0 ) then 
3730   f_riv_loc_c = f_riv_c(ji,jj) * friver_dep(jk,jmbathy) / fthk 
3731   fc_prod = fc_prod + f_riv_loc_c 
3732   endif 
3733   !! 
3734   !! benthic remineralisation 
3735   if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then 
3736   fc_prod = fc_prod + (f_benout_c(ji,jj) / fthk) 
3737   endif 
3738   if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then 
3739   fc_prod = fc_prod + (f_benout_ca(ji,jj) / fthk) 
3740   endif 
3741   !! 
3742   !! community respiration (does not include CaCO3 terms  obviously!) 
3743   fcomm_resp(ji,jj) = fcomm_resp(ji,jj) + fc_prod 
3744   !! 
3745   !! CaCO3 
3746   fc_prod = fc_prod  ftempca + freminca 
3747   !! 
3748   !! riverine flux 
3749   if ( jk .eq. 1 .and. jriver_c .gt. 0 ) then 
3750   fc_prod = fc_prod + f_riv_c(ji,jj) 
3751   endif 
3752   !! 
3753   btra(jpdic) = b0 * ( & 
3754   fc_prod + fc_cons ) 
3755   !! 
3756   fcar_cons(ji,jj) = fcar_cons(ji,jj) + ( fthk * ( & ! consumption of dissolved carbon 
3757   fc_cons ) ) 
3758   fcar_prod(ji,jj) = fcar_prod(ji,jj) + ( fthk * ( & ! production of dissolved carbon 
3759   fc_prod ) ) 
3760   !! 
3761   !! 
3762   !! alkalinity 
3763   fa_prod = 0.0 & 
3764   + (2.0 * freminca) ! CaCO3 dissolution 
3765   fa_cons = 0.0 & 
3766    (2.0 * ftempca) ! CaCO3 production 
3767   !! 
3768   !! riverine flux 
3769   if ( jriver_alk .gt. 0 ) then 
3770   f_riv_loc_alk = f_riv_alk(ji,jj) * friver_dep(jk,jmbathy) / fthk 
3771   fa_prod = fa_prod + f_riv_loc_alk 
3772   endif 
3773   !! 
3774   !! benthic remineralisation 
3775   if (jk.eq.jmbathy .and. jinorgben.eq.1 .and. ibenthic.eq.1) then 
3776   fa_prod = fa_prod + (2.0 * f_benout_ca(ji,jj) / fthk) 
3777   endif 
3778   !! 
3779   btra(jpalk) = b0 * ( & 
3780   fa_prod + fa_cons ) 
3781   !! 
3782   !! 
3783   !! oxygen (has protection at low O2 concentrations; OCMIP2 style) 
3784   fo2_prod = 0.0 & 
3785   + (xthetanit * fprn * zphn) & ! Pn primary production, N 
3786   + (xthetanit * fprd * zphd) & ! Pd primary production, N 
3787   + (xthetarem * xthetapn * fprn * zphn) & ! Pn primary production, C 
3788   + (xthetarem * xthetapd * fprd * zphd) ! Pd primary production, C 
3789   fo2_ncons = 0.0 & 
3790    (xthetanit * xphi * fgmipn) & ! Pn messy feeding remin., N 
3791    (xthetanit * xphi * fgmid) & ! D messy feeding remin., N 
3792    (xthetanit * xphi * fgmepn) & ! Pn messy feeding remin., N 
3793    (xthetanit * xphi * fgmepd) & ! Pd messy feeding remin., N 
3794    (xthetanit * xphi * fgmezmi) & ! Zi messy feeding remin., N 
3795    (xthetanit * xphi * fgmed) & ! D messy feeding remin., N 
3796    (xthetanit * fmiexcr) & ! microzoo excretion, N 
3797    (xthetanit * fmeexcr) & ! mesozoo excretion, N 
3798    (xthetanit * fdd) & ! slow detritus remin., N 
3799    (xthetanit * freminn) & ! fast detritus remin., N 
3800    (xthetanit * fdpn2) & ! Pn losses, N 
3801    (xthetanit * fdpd2) & ! Pd losses, N 
3802    (xthetanit * fdzmi2) & ! Zmi losses, N 
3803    (xthetanit * fdzme2) ! Zme losses, N 
3804   !! 
3805   !! benthic remineralisation 
3806   if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then 
3807   fo2_ncons = fo2_ncons  (xthetanit * f_benout_n(ji,jj) / fthk) 
3808   endif 
3809   fo2_ccons = 0.0 & 
3810    (xthetarem * xthetapn * xphi * fgmipn) & ! Pn messy feeding remin., C 
3811    (xthetarem * xphi * fgmidc) & ! D messy feeding remin., C 
3812    (xthetarem * xthetapn * xphi * fgmepn) & ! Pn messy feeding remin., C 
3813    (xthetarem * xthetapd * xphi * fgmepd) & ! Pd messy feeding remin., C 
3814    (xthetarem * xthetazmi * xphi * fgmezmi) & ! Zi messy feeding remin., C 
3815    (xthetarem * xphi * fgmedc) & ! D messy feeding remin., C 
3816    (xthetarem * fmiresp) & ! microzoo respiration, C 
3817    (xthetarem * fmeresp) & ! mesozoo respiration, C 
3818    (xthetarem * fddc) & ! slow detritus remin., C 
3819    (xthetarem * freminc) & ! fast detritus remin., C 
3820    (xthetarem * xthetapn * fdpn2) & ! Pn losses, C 
3821    (xthetarem * xthetapd * fdpd2) & ! Pd losses, C 
3822    (xthetarem * xthetazmi * fdzmi2) & ! Zmi losses, C 
3823    (xthetarem * xthetazme * fdzme2) ! Zme losses, C 
3824   !! 
3825   !! benthic remineralisation 
3826   if (jk.eq.jmbathy .and. jorgben.eq.1 .and. ibenthic.eq.1) then 
3827   fo2_ccons = fo2_ccons  (xthetarem * f_benout_c(ji,jj) / fthk) 
3828   endif 
3829   fo2_cons = fo2_ncons + fo2_ccons 
3830   !! 
3831   !! is this a suboxic zone? 
3832   if (zoxy.lt.xo2min) then ! deficient O2; production fluxes only 
3833   btra(jpoxy) = b0 * ( & 
3834   fo2_prod ) 
3835   foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod ) 
3836   foxy_anox(ji,jj) = foxy_anox(ji,jj) + ( fthk * fo2_cons ) 
3837   else ! sufficient O2; production + consumption fluxes 
3838   btra(jpoxy) = b0 * ( & 
3839   fo2_prod + fo2_cons ) 
3840   foxy_prod(ji,jj) = foxy_prod(ji,jj) + ( fthk * fo2_prod ) 
3841   foxy_cons(ji,jj) = foxy_cons(ji,jj) + ( fthk * fo2_cons ) 
3842   endif 
3843   !! 
3844   !! airsea fluxes (if this is the surface box) 
3845   if (jk.eq.1) then 
3846   !! 
3847   !! CO2 flux 
3848   btra(jpdic) = btra(jpdic) + (b0 * f_co2flux) 
3849   !! 
3850   !! O2 flux (mol/m3/s > mmol/m3/d) 
3851   btra(jpoxy) = btra(jpoxy) + (b0 * f_o2flux) 
3852   endif 
3853   # endif 
3854   
3855   # if defined key_debug_medusa 
3856   !! report state variable fluxes (not including fastsinking detritus) 
3857   if (idf.eq.1.AND.idfval.eq.1) then 
3858   IF (lwp) write (numout,*) '' 
3859   IF (lwp) write (numout,*) 'btra(jpchn)(',jk,') = ', btra(jpchn) 
3860   IF (lwp) write (numout,*) 'btra(jpchd)(',jk,') = ', btra(jpchd) 
3861   IF (lwp) write (numout,*) 'btra(jpphn)(',jk,') = ', btra(jpphn) 
3862   IF (lwp) write (numout,*) 'btra(jpphd)(',jk,') = ', btra(jpphd) 
3863   IF (lwp) write (numout,*) 'btra(jppds)(',jk,') = ', btra(jppds) 
3864   IF (lwp) write (numout,*) 'btra(jpzmi)(',jk,') = ', btra(jpzmi) 
3865   IF (lwp) write (numout,*) 'btra(jpzme)(',jk,') = ', btra(jpzme) 
3866   IF (lwp) write (numout,*) 'btra(jpdet)(',jk,') = ', btra(jpdet) 
3867   IF (lwp) write (numout,*) 'btra(jpdin)(',jk,') = ', btra(jpdin) 
3868   IF (lwp) write (numout,*) 'btra(jpsil)(',jk,') = ', btra(jpsil) 
3869   IF (lwp) write (numout,*) 'btra(jpfer)(',jk,') = ', btra(jpfer) 
3870   # if defined key_roam 
3871   IF (lwp) write (numout,*) 'btra(jpdtc)(',jk,') = ', btra(jpdtc) 
3872   IF (lwp) write (numout,*) 'btra(jpdic)(',jk,') = ', btra(jpdic) 
3873   IF (lwp) write (numout,*) 'btra(jpalk)(',jk,') = ', btra(jpalk) 
3874   IF (lwp) write (numout,*) 'btra(jpoxy)(',jk,') = ', btra(jpoxy) 
3875   # endif 
3876   endif 
3877   # endif 
3878   
3879   !! 
3880   !! Integrate calculated fluxes for mass balance 
3881   !! 
3882   !! 
3883   !! === nitrogen === 
3884   fflx_n(ji,jj) = fflx_n(ji,jj) + & 
3885   fthk * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) ) 
3886   !! === silicon === 
3887   fflx_si(ji,jj) = fflx_si(ji,jj) + & 
3888   fthk * ( btra(jppds) + btra(jpsil) ) 
3889   !! === iron === 
3890   fflx_fe(ji,jj) = fflx_fe(ji,jj) + & 
3891   fthk * ( ( xrfn * ( btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet)) ) + btra(jpfer) ) 
3892   # if defined key_roam 
3893   !! === carbon === 
3894   fflx_c(ji,jj) = fflx_c(ji,jj) + & 
3895   fthk * ( (xthetapn * btra(jpphn)) + (xthetapd * btra(jpphd)) + & 
3896   (xthetazmi * btra(jpzmi)) + (xthetazme * btra(jpzme)) + btra(jpdtc) + btra(jpdic) ) 
3897   !! === alkalinity === 
3898   fflx_a(ji,jj) = fflx_a(ji,jj) + & 
3899   fthk * ( btra(jpalk) ) 
3900   !! === oxygen === 
3901   fflx_o2(ji,jj) = fflx_o2(ji,jj) + & 
3902   fthk * ( btra(jpoxy) ) 
3903   # endif 
3904   
3905   !! 
3906   !! Apply calculated tracer fluxes 
3907   !! 
3908   !! 
3909   !! units: [unit of tracer] per second (fluxes are calculated above per day) 
3910   !! 
3911   ibio_switch = 1 
3912   # if defined key_gulf_finland 
3913   !! AXY (17/05/13): fudge in a Gulf of Finland correction; uses longitude 
3914   !! latitude range to establish if this is a Gulf of Finland 
3915   !! grid cell; if so, then BGC fluxes are ignored (though 
3916   !! still calculated); for reference, this is meant to be a 
3917   !! temporary fix to see if all of my problems can be done 
3918   !! away with if I switch off BGC fluxes in the Gulf of 
3919   !! Finland, which currently appears the source of trouble 
3920   if ( glamt(ji,jj).gt.24.7 .and. glamt(ji,jj).lt.27.8 .and. & 
3921   & gphit(ji,jj).gt.59.2 .and. gphit(ji,jj).lt.60.2 ) then 
3922   ibio_switch = 0 
3923   endif 
3924   # endif 
3925   if (ibio_switch.eq.1) then 
3926   tra(ji,jj,jk,jpchn) = tra(ji,jj,jk,jpchn) + (btra(jpchn) / 86400.) 
3927   tra(ji,jj,jk,jpchd) = tra(ji,jj,jk,jpchd) + (btra(jpchd) / 86400.) 
3928   tra(ji,jj,jk,jpphn) = tra(ji,jj,jk,jpphn) + (btra(jpphn) / 86400.) 
3929   tra(ji,jj,jk,jpphd) = tra(ji,jj,jk,jpphd) + (btra(jpphd) / 86400.) 
3930   tra(ji,jj,jk,jppds) = tra(ji,jj,jk,jppds) + (btra(jppds) / 86400.) 
3931   tra(ji,jj,jk,jpzmi) = tra(ji,jj,jk,jpzmi) + (btra(jpzmi) / 86400.) 
3932   tra(ji,jj,jk,jpzme) = tra(ji,jj,jk,jpzme) + (btra(jpzme) / 86400.) 
3933   tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + (btra(jpdet) / 86400.) 
3934   tra(ji,jj,jk,jpdin) = tra(ji,jj,jk,jpdin) + (btra(jpdin) / 86400.) 
3935   tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + (btra(jpsil) / 86400.) 
3936   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + (btra(jpfer) / 86400.) 
3937   # if defined key_roam 
3938   tra(ji,jj,jk,jpdtc) = tra(ji,jj,jk,jpdtc) + (btra(jpdtc) / 86400.) 
3939   tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + (btra(jpdic) / 86400.) 
3940   tra(ji,jj,jk,jpalk) = tra(ji,jj,jk,jpalk) + (btra(jpalk) / 86400.) 
3941   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + (btra(jpoxy) / 86400.) 
3942   # endif 
3943   endif 
3944   
3945   !! AXY (18/11/16): CMIP6 diagnostics 
3946   IF( med_diag%FBDDTALK%dgsave ) THEN 
3947   fbddtalk(ji,jj) = fbddtalk(ji,jj) + (btra(jpalk) * fthk) 
3948   ENDIF 
3949   IF( med_diag%FBDDTDIC%dgsave ) THEN 
3950   fbddtdic(ji,jj) = fbddtdic(ji,jj) + (btra(jpdic) * fthk) 
3951   ENDIF 
3952   IF( med_diag%FBDDTDIFE%dgsave ) THEN 
3953   fbddtdife(ji,jj) = fbddtdife(ji,jj) + (btra(jpfer) * fthk) 
3954   ENDIF 
3955   IF( med_diag%FBDDTDIN%dgsave ) THEN 
3956   fbddtdin(ji,jj) = fbddtdin(ji,jj) + (btra(jpdin) * fthk) 
3957   ENDIF 
3958   IF( med_diag%FBDDTDISI%dgsave ) THEN 
3959   fbddtdisi(ji,jj) = fbddtdisi(ji,jj) + (btra(jpsil) * fthk) 
3960   ENDIF 
3961   !! 
3962   IF( med_diag%BDDTALK3%dgsave ) THEN 
3963   bddtalk3(ji,jj,jk) = btra(jpalk) 
3964   ENDIF 
3965   IF( med_diag%BDDTDIC3%dgsave ) THEN 
3966   bddtdic3(ji,jj,jk) = btra(jpdic) 
3967   ENDIF 
3968   IF( med_diag%BDDTDIFE3%dgsave ) THEN 
3969   bddtdife3(ji,jj,jk) = btra(jpfer) 
3970   ENDIF 
3971   IF( med_diag%BDDTDIN3%dgsave ) THEN 
3972   bddtdin3(ji,jj,jk) = btra(jpdin) 
3973   ENDIF 
3974   IF( med_diag%BDDTDISI3%dgsave ) THEN 
3975   bddtdisi3(ji,jj,jk) = btra(jpsil) 
3976   ENDIF 
3977   
3978   # if defined key_debug_medusa 
3979   IF (lwp) write (numout,*) '' 
3980   IF (lwp) write (numout,*) 'trc_bio_medusa: end all calculations' 
3981   IF (lwp) write (numout,*) 'trc_bio_medusa: now outputs' 
3982   CALL flush(numout) 
3983   # endif 
3984   
3985   # if defined key_axy_nancheck 
3986   !! 
3987   !! Check calculated tracer fluxes 
3988   !! 
3989   !! 
3990   DO jn = 1,jptra 
3991   fq0 = btra(jn) 
3992   !! AXY (30/01/14): "isnan" problem on HECTOR 
3993   !! if (fq0 /= fq0 ) then 
3994   if ( ieee_is_nan( fq0 ) ) then 
3995   !! there's a NaN here 
3996   if (lwp) write(numout,*) 'NAN detected in btra(', ji, ',', & 
3997   & jj, ',', jk, ',', jn, ') at time', kt 
3998   CALL ctl_stop( 'trcbio_medusa, NAN in btra field' ) 
3999   endif 
4000   ENDDO 
4001   DO jn = 1,jptra 
4002   fq0 = tra(ji,jj,jk,jn) 
4003   !! AXY (30/01/14): "isnan" problem on HECTOR 
4004   !! if (fq0 /= fq0 ) then 
4005   if ( ieee_is_nan( fq0 ) ) then 
4006   !! there's a NaN here 
4007   if (lwp) write(numout,*) 'NAN detected in tra(', ji, ',', & 
4008   & jj, ',', jk, ',', jn, ') at time', kt 
4009   CALL ctl_stop( 'trcbio_medusa, NAN in tra field' ) 
4010   endif 
4011   ENDDO 
4012   CALL flush(numout) 
4013   # endif 
4014   
4015   !! 
4016   !! Check model conservation 
4017   !! these terms merely sum up the tendency terms of the relevant 
4018   !! state variables, which should sum to zero; the iron cycle is 
4019   !! complicated by fluxes that add (aeolian deposition and seafloor 
4020   !! remineralisation) and remove (scavenging) dissolved iron from 
4021   !! the model (i.e. the sum of iron fluxes is unlikely to be zero) 
4022   !! 
4023   !! 
4024   !! fnit0 = btra(jpphn) + btra(jpphd) + btra(jpzmi) + btra(jpzme) + btra(jpdet) + btra(jpdin) ! + ftempn 
4025   !! fsil0 = btra(jppds) + btra(jpsil) ! + ftempsi 
4026   !! ffer0 = (xrfn * fnit0) + btra(jpfer) 
4027   # if defined key_roam 
4028   !! fcar0 = 0. 
4029   !! falk0 = 0. 
4030   !! foxy0 = 0. 
4031   # endif 
4032   !! 
4033   !! if (kt/240*240.eq.kt) then 
4034   !! if (ji.eq.2.and.jj.eq.2.and.jk.eq.1) then 
4035   !! IF (lwp) write (*,*) '*******!MEDUSA Conservation!*******',kt 
4036   # if defined key_roam 
4037   !! IF (lwp) write (*,*) fnit0,fsil0,ffer0,fcar0,falk0,foxy0 
4038   # else 
4039   !! IF (lwp) write (*,*) fnit0,fsil0,ffer0 
4040   # endif 
4041   !! endif 
4042   !! endif 
4043   
4044   # if defined key_trc_diabio 
4045   !!====================================================================== 
4046   !! LOCAL GRID CELL DIAGNOSTICS 
4047   !!====================================================================== 
4048   !! 
4049   !! 
4050   !! Full diagnostics key_trc_diabio: 
4051   !! LOBSTER and PISCES support full diagnistics option key_trc_diabio 
4052   !! which gives an option of FULL output of biological sourses and sinks. 
4053   !! I cannot see any reason for doing this. If needed, it can be done 
4054   !! as shown below. 
4055   !! 
4056   !! 
4057   IF(lwp) WRITE(numout,*) ' MEDUSA does not support key_trc_diabio' 
4058   !! trbio(ji,jj,jk, 1) = fprn 
4059   # endif 
4060   
4061   IF( lk_iomput .AND. .NOT. ln_diatrc ) THEN 
4062   !! 
4063   !! Add in XML diagnostics stuff 
4064   !! 
4065   !! 
4066   !! ** 2D diagnostics 
4067   # if defined key_debug_medusa 
4068   IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ijjjjk loop' 
4069   CALL flush(numout) 
4070   # endif 
4071   IF ( med_diag%PRN%dgsave ) THEN 
4072   fprn2d(ji,jj) = fprn2d(ji,jj) + (fprn * zphn * fthk) 
4073   ENDIF 
4074   IF ( med_diag%MPN%dgsave ) THEN 
4075   fdpn2d(ji,jj) = fdpn2d(ji,jj) + (fdpn * fthk) 
4076   ENDIF 
4077   IF ( med_diag%PRD%dgsave ) THEN 
4078   fprd2d(ji,jj) = fprd2d(ji,jj) + (fprd * zphd * fthk) 
4079   ENDIF 
4080   IF( med_diag%MPD%dgsave ) THEN 
4081   fdpd2d(ji,jj) = fdpd2d(ji,jj) + (fdpd * fthk) 
4082   ENDIF 
4083   ! IF( med_diag%DSED%dgsave ) THEN 
4084   ! CALL iom_put( "DSED" , ftot_n ) 
4085   ! ENDIF 
4086   IF( med_diag%OPAL%dgsave ) THEN 
4087   fprds2d(ji,jj) = fprds2d(ji,jj) + (fprds * zpds * fthk) 
4088   ENDIF 
4089   IF( med_diag%OPALDISS%dgsave ) THEN 
4090   fsdiss2d(ji,jj) = fsdiss2d(ji,jj) + (fsdiss * fthk) 
4091   ENDIF 
4092   IF( med_diag%GMIPn%dgsave ) THEN 
4093   fgmipn2d(ji,jj) = fgmipn2d(ji,jj) + (fgmipn * fthk) 
4094   ENDIF 
4095   IF( med_diag%GMID%dgsave ) THEN 
4096   fgmid2d(ji,jj) = fgmid2d(ji,jj) + (fgmid * fthk) 
4097   ENDIF 
4098   IF( med_diag%MZMI%dgsave ) THEN 
4099   fdzmi2d(ji,jj) = fdzmi2d(ji,jj) + (fdzmi * fthk) 
4100   ENDIF 
4101   IF( med_diag%GMEPN%dgsave ) THEN 
4102   fgmepn2d(ji,jj) = fgmepn2d(ji,jj) + (fgmepn * fthk) 
4103   ENDIF 
4104   IF( med_diag%GMEPD%dgsave ) THEN 
4105   fgmepd2d(ji,jj) = fgmepd2d(ji,jj) + (fgmepd * fthk) 
4106   ENDIF 
4107   IF( med_diag%GMEZMI%dgsave ) THEN 
4108   fgmezmi2d(ji,jj) = fgmezmi2d(ji,jj) + (fgmezmi * fthk) 
4109   ENDIF 
4110   IF( med_diag%GMED%dgsave ) THEN 
4111   fgmed2d(ji,jj) = fgmed2d(ji,jj) + (fgmed * fthk) 
4112   ENDIF 
4113   IF( med_diag%MZME%dgsave ) THEN 
4114   fdzme2d(ji,jj) = fdzme2d(ji,jj) + (fdzme * fthk) 
4115   ENDIF 
4116   ! IF( med_diag%DEXP%dgsave ) THEN 
4117   ! CALL iom_put( "DEXP" , ftot_n ) 
4118   ! ENDIF 
4119   IF( med_diag%DETN%dgsave ) THEN 
4120   fslown2d(ji,jj) = fslown2d(ji,jj) + (fslown * fthk) 
4121   ENDIF 
4122   IF( med_diag%MDET%dgsave ) THEN 
4123   fdd2d(ji,jj) = fdd2d(ji,jj) + (fdd * fthk) 
4124   ENDIF 
4125   IF( med_diag%AEOLIAN%dgsave ) THEN 
4126   ffetop2d(ji,jj) = ffetop2d(ji,jj) + (ffetop * fthk) 
4127   ENDIF 
4128   IF( med_diag%BENTHIC%dgsave ) THEN 
4129   ffebot2d(ji,jj) = ffebot2d(ji,jj) + (ffebot * fthk) 
4130   ENDIF 
4131   IF( med_diag%SCAVENGE%dgsave ) THEN 
4132   ffescav2d(ji,jj) = ffescav2d(ji,jj) + (ffescav * fthk) 
4133   ENDIF 
4134   IF( med_diag%PN_JLIM%dgsave ) THEN 
4135   ! fjln2d(ji,jj) = fjln2d(ji,jj) + (fjln * zphn * fthk) 
4136   fjln2d(ji,jj) = fjln2d(ji,jj) + (fjlim_pn * zphn * fthk) 
4137   ENDIF 
4138   IF( med_diag%PN_NLIM%dgsave ) THEN 
4139   fnln2d(ji,jj) = fnln2d(ji,jj) + (fnln * zphn * fthk) 
4140   ENDIF 
4141   IF( med_diag%PN_FELIM%dgsave ) THEN 
4142   ffln2d(ji,jj) = ffln2d(ji,jj) + (ffln * zphn * fthk) 
4143   ENDIF 
4144   IF( med_diag%PD_JLIM%dgsave ) THEN 
4145   ! fjld2d(ji,jj) = fjld2d(ji,jj) + (fjld * zphd * fthk) 
4146   fjld2d(ji,jj) = fjld2d(ji,jj) + (fjlim_pd * zphd * fthk) 
4147   ENDIF 
4148   IF( med_diag%PD_NLIM%dgsave ) THEN 
4149   fnld2d(ji,jj) = fnld2d(ji,jj) + (fnld * zphd * fthk) 
4150   ENDIF 
4151   IF( med_diag%PD_FELIM%dgsave ) THEN 
4152   ffld2d(ji,jj) = ffld2d(ji,jj) + (ffld * zphd * fthk) 
4153   ENDIF 
4154   IF( med_diag%PD_SILIM%dgsave ) THEN 
4155   fsld2d2(ji,jj) = fsld2d2(ji,jj) + (fsld2 * zphd * fthk) 
4156   ENDIF 
4157   IF( med_diag%PDSILIM2%dgsave ) THEN 
4158   fsld2d(ji,jj) = fsld2d(ji,jj) + (fsld * zphd * fthk) 
4159   ENDIF 
4160   !! 
4161   IF( med_diag%TOTREG_N%dgsave ) THEN 
4162   fregen2d(ji,jj) = fregen2d(ji,jj) + fregen 
4163   ENDIF 
4164   IF( med_diag%TOTRG_SI%dgsave ) THEN 
4165   fregensi2d(ji,jj) = fregensi2d(ji,jj) + fregensi 
4166   ENDIF 
4167   !! 
4168   IF( med_diag%FASTN%dgsave ) THEN 
4169   ftempn2d(ji,jj) = ftempn2d(ji,jj) + (ftempn * fthk) 
4170   ENDIF 
4171   IF( med_diag%FASTSI%dgsave ) THEN 
4172   ftempsi2d(ji,jj) = ftempsi2d(ji,jj) + (ftempsi * fthk) 
4173   ENDIF 
4174   IF( med_diag%FASTFE%dgsave ) THEN 
4175   ftempfe2d(ji,jj) =ftempfe2d(ji,jj) + (ftempfe * fthk) 
4176   ENDIF 
4177   IF( med_diag%FASTC%dgsave ) THEN 
4178   ftempc2d(ji,jj) = ftempc2d(ji,jj) + (ftempc * fthk) 
4179   ENDIF 
4180   IF( med_diag%FASTCA%dgsave ) THEN 
4181   ftempca2d(ji,jj) = ftempca2d(ji,jj) + (ftempca * fthk) 
4182   ENDIF 
4183   !! 
4184   IF( med_diag%REMINN%dgsave ) THEN 
4185   freminn2d(ji,jj) = freminn2d(ji,jj) + (freminn * fthk) 
4186   ENDIF 
4187   IF( med_diag%REMINSI%dgsave ) THEN 
4188   freminsi2d(ji,jj) = freminsi2d(ji,jj) + (freminsi * fthk) 
4189   ENDIF 
4190   IF( med_diag%REMINFE%dgsave ) THEN 
4191   freminfe2d(ji,jj)= freminfe2d(ji,jj) + (freminfe * fthk) 
4192   ENDIF 
4193   IF( med_diag%REMINC%dgsave ) THEN 
4194   freminc2d(ji,jj) = freminc2d(ji,jj) + (freminc * fthk) 
4195   ENDIF 
4196   IF( med_diag%REMINCA%dgsave ) THEN 
4197   freminca2d(ji,jj) = freminca2d(ji,jj) + (freminca * fthk) 
4198   ENDIF 
4199   !! 
4200   # if defined key_roam 
4201   !! 
4202   !! AXY (09/11/16): CMIP6 diagnostics 
4203   IF( med_diag%FD_NIT3%dgsave ) THEN 
4204   fd_nit3(ji,jj,jk) = ffastn(ji,jj) 
4205   ENDIF 
4206   IF( med_diag%FD_SIL3%dgsave ) THEN 
4207   fd_sil3(ji,jj,jk) = ffastsi(ji,jj) 
4208   ENDIF 
4209   IF( med_diag%FD_CAR3%dgsave ) THEN 
4210   fd_car3(ji,jj,jk) = ffastc(ji,jj) 
4211   ENDIF 
4212   IF( med_diag%FD_CAL3%dgsave ) THEN 
4213   fd_cal3(ji,jj,jk) = ffastca(ji,jj) 
4214   ENDIF 
4215   !! 
4216   IF (jk.eq.i0100) THEN 
4217   IF( med_diag%RR_0100%dgsave ) THEN 
4218   ffastca2d(ji,jj) = & 
4219   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) 
4220   ENDIF 
4221   ELSE IF (jk.eq.i0500) THEN 
4222   IF( med_diag%RR_0500%dgsave ) THEN 
4223   ffastca2d(ji,jj) = & 
4224   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) 
4225   ENDIF 
4226   ELSE IF (jk.eq.i1000) THEN 
4227   IF( med_diag%RR_1000%dgsave ) THEN 
4228   ffastca2d(ji,jj) = & 
4229   ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) 
4230   ENDIF 
4231   ELSE IF (jk.eq.jmbathy) THEN 
4232   IF( med_diag%IBEN_N%dgsave ) THEN 
4233   iben_n2d(ji,jj) = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) 
4234   ENDIF 
4235   IF( med_diag%IBEN_FE%dgsave ) THEN 
4236   iben_fe2d(ji,jj) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj) 
4237   ENDIF 
4238   IF( med_diag%IBEN_C%dgsave ) THEN 
4239   iben_c2d(ji,jj) = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) 
4240   ENDIF 
4241   IF( med_diag%IBEN_SI%dgsave ) THEN 
4242   iben_si2d(ji,jj) = f_fbenin_si(ji,jj) 
4243   ENDIF 
4244   IF( med_diag%IBEN_CA%dgsave ) THEN 
4245   iben_ca2d(ji,jj) = f_fbenin_ca(ji,jj) 
4246   ENDIF 
4247   IF( med_diag%OBEN_N%dgsave ) THEN 
4248   oben_n2d(ji,jj) = f_benout_n(ji,jj) 
4249   ENDIF 
4250   IF( med_diag%OBEN_FE%dgsave ) THEN 
4251   oben_fe2d(ji,jj) = f_benout_fe(ji,jj) 
4252   ENDIF 
4253   IF( med_diag%OBEN_C%dgsave ) THEN 
4254   oben_c2d(ji,jj) = f_benout_c(ji,jj) 
4255   ENDIF 
4256   IF( med_diag%OBEN_SI%dgsave ) THEN 
4257   oben_si2d(ji,jj) = f_benout_si(ji,jj) 
4258   ENDIF 
4259   IF( med_diag%OBEN_CA%dgsave ) THEN 
4260   oben_ca2d(ji,jj) = f_benout_ca(ji,jj) 
4261   ENDIF 
4262   IF( med_diag%SFR_OCAL%dgsave ) THEN 
4263   sfr_ocal2d(ji,jj) = f3_omcal(ji,jj,jk) 
4264   ENDIF 
4265   IF( med_diag%SFR_OARG%dgsave ) THEN 
4266   sfr_oarg2d(ji,jj) = f3_omarg(ji,jj,jk) 
4267   ENDIF 
4268   IF( med_diag%LYSO_CA%dgsave ) THEN 
4269   lyso_ca2d(ji,jj) = f_benout_lyso_ca(ji,jj) 
4270   ENDIF 
4271   ENDIF 
4272   !! end bathy1 diags 
4273   !! 
4274   IF( med_diag%RIV_N%dgsave ) THEN 
4275   rivn2d(ji,jj) = rivn2d(ji,jj) + (f_riv_loc_n * fthk) 
4276   ENDIF 
4277   IF( med_diag%RIV_SI%dgsave ) THEN 
4278   rivsi2d(ji,jj) = rivsi2d(ji,jj) + (f_riv_loc_si * fthk) 
4279   ENDIF 
4280   IF( med_diag%RIV_C%dgsave ) THEN 
4281   rivc2d(ji,jj) = rivc2d(ji,jj) + (f_riv_loc_c * fthk) 
4282   ENDIF 
4283   IF( med_diag%RIV_ALK%dgsave ) THEN 
4284   rivalk2d(ji,jj) = rivalk2d(ji,jj) + (f_riv_loc_alk * fthk) 
4285   ENDIF 
4286   IF( med_diag%DETC%dgsave ) THEN 
4287   fslowc2d(ji,jj) = fslowc2d(ji,jj) + (fslowc * fthk) 
4288   ENDIF 
4289   !! 
4290   !! 
4291   !! 
4292   IF( med_diag%PN_LLOSS%dgsave ) THEN 
4293   fdpn22d(ji,jj) = fdpn22d(ji,jj) + (fdpn2 * fthk) 
4294   ENDIF 
4295   IF( med_diag%PD_LLOSS%dgsave ) THEN 
4296   fdpd22d(ji,jj) = fdpd22d(ji,jj) + (fdpd2 * fthk) 
4297   ENDIF 
4298   IF( med_diag%ZI_LLOSS%dgsave ) THEN 
4299   fdzmi22d(ji,jj) = fdzmi22d(ji,jj) + (fdzmi2 * fthk) 
4300   ENDIF 
4301   IF( med_diag%ZE_LLOSS%dgsave ) THEN 
4302   fdzme22d(ji,jj) = fdzme22d(ji,jj) + (fdzme2 * fthk) 
4303   ENDIF 
4304   IF( med_diag%ZI_MES_N%dgsave ) THEN 
4305   zimesn2d(ji,jj) = zimesn2d(ji,jj) + & 
4306   (xphi * (fgmipn + fgmid) * fthk) 
4307   ENDIF 
4308   IF( med_diag%ZI_MES_D%dgsave ) THEN 
4309   zimesd2d(ji,jj) = zimesd2d(ji,jj) + & 
4310   ((1.  xbetan) * finmi * fthk) 
4311   ENDIF 
4312   IF( med_diag%ZI_MES_C%dgsave ) THEN 
4313   zimesc2d(ji,jj) = zimesc2d(ji,jj) + & 
4314   (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk) 
4315   ENDIF 
4316   IF( med_diag%ZI_MESDC%dgsave ) THEN 
4317   zimesdc2d(ji,jj) = zimesdc2d(ji,jj) + & 
4318   ((1.  xbetac) * ficmi * fthk) 
4319   ENDIF 
4320   IF( med_diag%ZI_EXCR%dgsave ) THEN 
4321   ziexcr2d(ji,jj) = ziexcr2d(ji,jj) + (fmiexcr * fthk) 
4322   ENDIF 
4323   IF( med_diag%ZI_RESP%dgsave ) THEN 
4324   ziresp2d(ji,jj) = ziresp2d(ji,jj) + (fmiresp * fthk) 
4325   ENDIF 
4326   IF( med_diag%ZI_GROW%dgsave ) THEN 
4327   zigrow2d(ji,jj) = zigrow2d(ji,jj) + (fmigrow * fthk) 
4328   ENDIF 
4329   IF( med_diag%ZE_MES_N%dgsave ) THEN 
4330   zemesn2d(ji,jj) = zemesn2d(ji,jj) + & 
4331   (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk) 
4332   ENDIF 
4333   IF( med_diag%ZE_MES_D%dgsave ) THEN 
4334   zemesd2d(ji,jj) = zemesd2d(ji,jj) + & 
4335   ((1.  xbetan) * finme * fthk) 
4336   ENDIF 
4337   IF( med_diag%ZE_MES_C%dgsave ) THEN 
4338   zemesc2d(ji,jj) = zemesc2d(ji,jj) + & 
4339   (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & 
4340   (xthetazmi * fgmezmi) + fgmedc) * fthk) 
4341   ENDIF 
4342   IF( med_diag%ZE_MESDC%dgsave ) THEN 
4343   zemesdc2d(ji,jj) = zemesdc2d(ji,jj) + & 
4344   ((1.  xbetac) * ficme * fthk) 
4345   ENDIF 
4346   IF( med_diag%ZE_EXCR%dgsave ) THEN 
4347   zeexcr2d(ji,jj) = zeexcr2d(ji,jj) + (fmeexcr * fthk) 
4348   ENDIF 
4349   IF( med_diag%ZE_RESP%dgsave ) THEN 
4350   zeresp2d(ji,jj) = zeresp2d(ji,jj) + (fmeresp * fthk) 
4351   ENDIF 
4352   IF( med_diag%ZE_GROW%dgsave ) THEN 
4353   zegrow2d(ji,jj) = zegrow2d(ji,jj) + (fmegrow * fthk) 
4354   ENDIF 
4355   IF( med_diag%MDETC%dgsave ) THEN 
4356   mdetc2d(ji,jj) = mdetc2d(ji,jj) + (fddc * fthk) 
4357   ENDIF 
4358   IF( med_diag%GMIDC%dgsave ) THEN 
4359   gmidc2d(ji,jj) = gmidc2d(ji,jj) + (fgmidc * fthk) 
4360   ENDIF 
4361   IF( med_diag%GMEDC%dgsave ) THEN 
4362   gmedc2d(ji,jj) = gmedc2d(ji,jj) + (fgmedc * fthk) 
4363   ENDIF 
4364   !! 
4365   # endif 
4366   !! 
4367   !! ** 3D diagnostics 
4368   IF( med_diag%TPP3%dgsave ) THEN 
4369   tpp3d(ji,jj,jk) = (fprn * zphn) + (fprd * zphd) 
4370   !CALL iom_put( "TPP3" , tpp3d ) 
4371   ENDIF 
4372   IF( med_diag%TPPD3%dgsave ) THEN 
4373   tppd3(ji,jj,jk) = (fprd * zphd) 
4374   ENDIF 
4375   
4376   IF( med_diag%REMIN3N%dgsave ) THEN 
4377   remin3dn(ji,jj,jk) = fregen + (freminn * fthk) !! remineralisation 
4378   !CALL iom_put( "REMIN3N" , remin3dn ) 
4379   ENDIF 
4380   !! IF( med_diag%PH3%dgsave ) THEN 
4381   !! CALL iom_put( "PH3" , f3_pH ) 
4382   !! ENDIF 
4383   !! IF( med_diag%OM_CAL3%dgsave ) THEN 
4384   !! CALL iom_put( "OM_CAL3" , f3_omcal ) 
4385   !! ENDIF 
4386   !! 
4387   !! AXY (09/11/16): CMIP6 diagnostics 
4388   IF ( med_diag%DCALC3%dgsave ) THEN 
4389   dcalc3(ji,jj,jk) = freminca 
4390   ENDIF 
4391   IF ( med_diag%FEDISS3%dgsave ) THEN 
4392   fediss3(ji,jj,jk) = ffetop 
4393   ENDIF 
4394   IF ( med_diag%FESCAV3%dgsave ) THEN 
4395   fescav3(ji,jj,jk) = ffescav 
4396   ENDIF 
4397   IF ( med_diag%MIGRAZP3%dgsave ) THEN 
4398   migrazp3(ji,jj,jk) = fgmipn * xthetapn 
4399   ENDIF 
4400   IF ( med_diag%MIGRAZD3%dgsave ) THEN 
4401   migrazd3(ji,jj,jk) = fgmidc 
4402   ENDIF 
4403   IF ( med_diag%MEGRAZP3%dgsave ) THEN 
4404   megrazp3(ji,jj,jk) = (fgmepn * xthetapn) + (fgmepd * xthetapd) 
4405   ENDIF 
4406   IF ( med_diag%MEGRAZD3%dgsave ) THEN 
4407   megrazd3(ji,jj,jk) = fgmedc 
4408   ENDIF 
4409   IF ( med_diag%MEGRAZZ3%dgsave ) THEN 
4410   megrazz3(ji,jj,jk) = (fgmezmi * xthetazmi) 
4411   ENDIF 
4412   IF ( med_diag%PBSI3%dgsave ) THEN 
4413   pbsi3(ji,jj,jk) = (fprds * zpds) 
4414   ENDIF 
4415   IF ( med_diag%PCAL3%dgsave ) THEN 
4416   pcal3(ji,jj,jk) = ftempca 
4417   ENDIF 
4418   IF ( med_diag%REMOC3%dgsave ) THEN 
4419   remoc3(ji,jj,jk) = freminc 
4420   ENDIF 
4421   IF ( med_diag%PNLIMJ3%dgsave ) THEN 
4422   ! pnlimj3(ji,jj,jk) = fjln 
4423   pnlimj3(ji,jj,jk) = fjlim_pn 
4424   ENDIF 
4425   IF ( med_diag%PNLIMN3%dgsave ) THEN 
4426   pnlimn3(ji,jj,jk) = fnln 
4427   ENDIF 
4428   IF ( med_diag%PNLIMFE3%dgsave ) THEN 
4429   pnlimfe3(ji,jj,jk) = ffln 
4430   ENDIF 
4431   IF ( med_diag%PDLIMJ3%dgsave ) THEN 
4432   ! pdlimj3(ji,jj,jk) = fjld 
4433   pdlimj3(ji,jj,jk) = fjlim_pd 
4434   ENDIF 
4435   IF ( med_diag%PDLIMN3%dgsave ) THEN 
4436   pdlimn3(ji,jj,jk) = fnld 
4437   ENDIF 
4438   IF ( med_diag%PDLIMFE3%dgsave ) THEN 
4439   pdlimfe3(ji,jj,jk) = ffld 
4440   ENDIF 
4441   IF ( med_diag%PDLIMSI3%dgsave ) THEN 
4442   pdlimsi3(ji,jj,jk) = fsld2 
4443   ENDIF 
4444   !! 
4445   !! ** Without using iom_use 
4446   ELSE IF( ln_diatrc ) THEN 
4447   # if defined key_debug_medusa 
4448   IF (lwp) write (numout,*) 'trc_bio_medusa: diag in ijjjjk ln_diatrc' 
4449   CALL flush(numout) 
4450   # endif 
4451   !! 
4452   !! Prepare 2D diagnostics 
4453   !! 
4454   !! 
4455   !! if ((kt / 240*240).eq.kt) then 
4456   !! IF (lwp) write (*,*) '*******!MEDUSA DIAADD!*******',kt 
4457   !! endif 
4458   trc2d(ji,jj,1) = ftot_n(ji,jj) !! nitrogen inventory 
4459   trc2d(ji,jj,2) = ftot_si(ji,jj) !! silicon inventory 
4460   trc2d(ji,jj,3) = ftot_fe(ji,jj) !! iron inventory 
4461   trc2d(ji,jj,4) = trc2d(ji,jj,4) + (fprn * zphn * fthk) !! nondiatom production 
4462   trc2d(ji,jj,5) = trc2d(ji,jj,5) + (fdpn * fthk) !! nondiatom nongrazing losses 
4463   trc2d(ji,jj,6) = trc2d(ji,jj,6) + (fprd * zphd * fthk) !! diatom production 
4464   trc2d(ji,jj,7) = trc2d(ji,jj,7) + (fdpd * fthk) !! diatom nongrazing losses 
4465   !! diagnostic field 8 is (ostensibly) supplied by trcsed.F 
4466   trc2d(ji,jj,9) = trc2d(ji,jj,9) + (fprds * zpds * fthk) !! diatom silicon production 
4467   trc2d(ji,jj,10) = trc2d(ji,jj,10) + (fsdiss * fthk) !! diatom silicon dissolution 
4468   trc2d(ji,jj,11) = trc2d(ji,jj,11) + (fgmipn * fthk) !! microzoo grazing on nondiatoms 
4469   trc2d(ji,jj,12) = trc2d(ji,jj,12) + (fgmid * fthk) !! microzoo grazing on detrital nitrogen 
4470   trc2d(ji,jj,13) = trc2d(ji,jj,13) + (fdzmi * fthk) !! microzoo nongrazing losses 
4471   trc2d(ji,jj,14) = trc2d(ji,jj,14) + (fgmepn * fthk) !! mesozoo grazing on nondiatoms 
4472   trc2d(ji,jj,15) = trc2d(ji,jj,15) + (fgmepd * fthk) !! mesozoo grazing on diatoms 
4473   trc2d(ji,jj,16) = trc2d(ji,jj,16) + (fgmezmi * fthk) !! mesozoo grazing on microzoo 
4474   trc2d(ji,jj,17) = trc2d(ji,jj,17) + (fgmed * fthk) !! mesozoo grazing on detrital nitrogen 
4475   trc2d(ji,jj,18) = trc2d(ji,jj,18) + (fdzme * fthk) !! mesozoo nongrazing losses 
4476   !! diagnostic field 19 is (ostensibly) supplied by trcexp.F 
4477   trc2d(ji,jj,20) = trc2d(ji,jj,20) + (fslown * fthk) !! slow sinking detritus N production 
4478   trc2d(ji,jj,21) = trc2d(ji,jj,21) + (fdd * fthk) !! detrital remineralisation 
4479   trc2d(ji,jj,22) = trc2d(ji,jj,22) + (ffetop * fthk) !! aeolian iron addition 
4480   trc2d(ji,jj,23) = trc2d(ji,jj,23) + (ffebot * fthk) !! seafloor iron addition 
4481   trc2d(ji,jj,24) = trc2d(ji,jj,24) + (ffescav * fthk) !! "free" iron scavenging 
4482   trc2d(ji,jj,25) = trc2d(ji,jj,25) + (fjlim_pn * zphn * fthk) !! nondiatom J limitation term 
4483   trc2d(ji,jj,26) = trc2d(ji,jj,26) + (fnln * zphn * fthk) !! nondiatom N limitation term 
4484   trc2d(ji,jj,27) = trc2d(ji,jj,27) + (ffln * zphn * fthk) !! nondiatom Fe limitation term 
4485   trc2d(ji,jj,28) = trc2d(ji,jj,28) + (fjlim_pd * zphd * fthk) !! diatom J limitation term 
4486   trc2d(ji,jj,29) = trc2d(ji,jj,29) + (fnld * zphd * fthk) !! diatom N limitation term 
4487   trc2d(ji,jj,30) = trc2d(ji,jj,30) + (ffld * zphd * fthk) !! diatom Fe limitation term 
4488   trc2d(ji,jj,31) = trc2d(ji,jj,31) + (fsld2 * zphd * fthk) !! diatom Si limitation term 
4489   trc2d(ji,jj,32) = trc2d(ji,jj,32) + (fsld * zphd * fthk) !! diatom Si uptake limitation term 
4490   if (jk.eq.i0100) trc2d(ji,jj,33) = fslownflux(ji,jj) !! slow detritus flux at 100 m 
4491   if (jk.eq.i0200) trc2d(ji,jj,34) = fslownflux(ji,jj) !! slow detritus flux at 200 m 
4492   if (jk.eq.i0500) trc2d(ji,jj,35) = fslownflux(ji,jj) !! slow detritus flux at 500 m 
4493   if (jk.eq.i1000) trc2d(ji,jj,36) = fslownflux(ji,jj) !! slow detritus flux at 1000 m 
4494   trc2d(ji,jj,37) = trc2d(ji,jj,37) + fregen !! nonfast N full column regeneration 
4495   trc2d(ji,jj,38) = trc2d(ji,jj,38) + fregensi !! nonfast Si full column regeneration 
4496   if (jk.eq.i0100) trc2d(ji,jj,39) = trc2d(ji,jj,37) !! nonfast N regeneration to 100 m 
4497   if (jk.eq.i0200) trc2d(ji,jj,40) = trc2d(ji,jj,37) !! nonfast N regeneration to 200 m 
4498   if (jk.eq.i0500) trc2d(ji,jj,41) = trc2d(ji,jj,37) !! nonfast N regeneration to 500 m 
4499   if (jk.eq.i1000) trc2d(ji,jj,42) = trc2d(ji,jj,37) !! nonfast N regeneration to 1000 m 
4500   trc2d(ji,jj,43) = trc2d(ji,jj,43) + (ftempn * fthk) !! fast sinking detritus N production 
4501   trc2d(ji,jj,44) = trc2d(ji,jj,44) + (ftempsi * fthk) !! fast sinking detritus Si production 
4502   trc2d(ji,jj,45) = trc2d(ji,jj,45) + (ftempfe * fthk) !! fast sinking detritus Fe production 
4503   trc2d(ji,jj,46) = trc2d(ji,jj,46) + (ftempc * fthk) !! fast sinking detritus C production 
4504   trc2d(ji,jj,47) = trc2d(ji,jj,47) + (ftempca * fthk) !! fast sinking detritus CaCO3 production 
4505   if (jk.eq.i0100) trc2d(ji,jj,48) = ffastn(ji,jj) !! fast detritus N flux at 100 m 
4506   if (jk.eq.i0200) trc2d(ji,jj,49) = ffastn(ji,jj) !! fast detritus N flux at 200 m 
4507   if (jk.eq.i0500) trc2d(ji,jj,50) = ffastn(ji,jj) !! fast detritus N flux at 500 m 
4508   if (jk.eq.i1000) trc2d(ji,jj,51) = ffastn(ji,jj) !! fast detritus N flux at 1000 m 
4509   if (jk.eq.i0100) trc2d(ji,jj,52) = fregenfast(ji,jj) !! N regeneration to 100 m 
4510   if (jk.eq.i0200) trc2d(ji,jj,53) = fregenfast(ji,jj) !! N regeneration to 200 m 
4511   if (jk.eq.i0500) trc2d(ji,jj,54) = fregenfast(ji,jj) !! N regeneration to 500 m 
4512   if (jk.eq.i1000) trc2d(ji,jj,55) = fregenfast(ji,jj) !! N regeneration to 1000 m 
4513   if (jk.eq.i0100) trc2d(ji,jj,56) = ffastsi(ji,jj) !! fast detritus Si flux at 100 m 
4514   if (jk.eq.i0200) trc2d(ji,jj,57) = ffastsi(ji,jj) !! fast detritus Si flux at 200 m 
4515   if (jk.eq.i0500) trc2d(ji,jj,58) = ffastsi(ji,jj) !! fast detritus Si flux at 500 m 
4516   if (jk.eq.i1000) trc2d(ji,jj,59) = ffastsi(ji,jj) !! fast detritus Si flux at 1000 m 
4517   if (jk.eq.i0100) trc2d(ji,jj,60) = fregenfastsi(ji,jj) !! Si regeneration to 100 m 
4518   if (jk.eq.i0200) trc2d(ji,jj,61) = fregenfastsi(ji,jj) !! Si regeneration to 200 m 
4519   if (jk.eq.i0500) trc2d(ji,jj,62) = fregenfastsi(ji,jj) !! Si regeneration to 500 m 
4520   if (jk.eq.i1000) trc2d(ji,jj,63) = fregenfastsi(ji,jj) !! Si regeneration to 1000 m 
4521   trc2d(ji,jj,64) = trc2d(ji,jj,64) + (freminn * fthk) !! sum of fastsinking N fluxes 
4522   trc2d(ji,jj,65) = trc2d(ji,jj,65) + (freminsi * fthk) !! sum of fastsinking Si fluxes 
4523   trc2d(ji,jj,66) = trc2d(ji,jj,66) + (freminfe * fthk) !! sum of fastsinking Fe fluxes 
4524   trc2d(ji,jj,67) = trc2d(ji,jj,67) + (freminc * fthk) !! sum of fastsinking C fluxes 
4525   trc2d(ji,jj,68) = trc2d(ji,jj,68) + (freminca * fthk) !! sum of fastsinking Ca fluxes 
4526   if (jk.eq.jmbathy) then 
4527   trc2d(ji,jj,69) = fsedn(ji,jj) !! N sedimentation flux 
4528   trc2d(ji,jj,70) = fsedsi(ji,jj) !! Si sedimentation flux 
4529   trc2d(ji,jj,71) = fsedfe(ji,jj) !! Fe sedimentation flux 
4530   trc2d(ji,jj,72) = fsedc(ji,jj) !! C sedimentation flux 
4531   trc2d(ji,jj,73) = fsedca(ji,jj) !! Ca sedimentation flux 
4532   endif 
4533   if (jk.eq.1) trc2d(ji,jj,74) = qsr(ji,jj) 
4534   if (jk.eq.1) trc2d(ji,jj,75) = xpar(ji,jj,jk) 
4535   !! if (jk.eq.1) trc2d(ji,jj,75) = real(iters) 
4536   !! diagnostic fields 76 to 80 calculated below 
4537   trc2d(ji,jj,81) = trc2d(ji,jj,81) + fprn_ml(ji,jj) !! mixed layer nondiatom production 
4538   trc2d(ji,jj,82) = trc2d(ji,jj,82) + fprd_ml(ji,jj) !! mixed layer diatom production 
4539   # if defined key_gulf_finland 
4540   if (jk.eq.1) trc2d(ji,jj,83) = real(ibio_switch) !! Gulf of Finland check 
4541   # else 
4542   trc2d(ji,jj,83) = ocal_ccd(ji,jj) !! calcite CCD depth 
4543   # endif 
4544   trc2d(ji,jj,84) = fccd(ji,jj) !! last model level above calcite CCD depth 
4545   if (jk.eq.1) trc2d(ji,jj,85) = xFree(ji,jj) !! surface "free" iron 
4546   if (jk.eq.i0200) trc2d(ji,jj,86) = xFree(ji,jj) !! "free" iron at 100 m 
4547   if (jk.eq.i0200) trc2d(ji,jj,87) = xFree(ji,jj) !! "free" iron at 200 m 
4548   if (jk.eq.i0500) trc2d(ji,jj,88) = xFree(ji,jj) !! "free" iron at 500 m 
4549   if (jk.eq.i1000) trc2d(ji,jj,89) = xFree(ji,jj) !! "free" iron at 1000 m 
4550   !! AXY (27/06/12): extract "euphotic depth" 
4551   if (jk.eq.1) trc2d(ji,jj,90) = xze(ji,jj) 
4552   !! 
4553   # if defined key_roam 
4554   !! ROAM provisionally has access to a further 20 2D diagnostics 
4555   if (jk .eq. 1) then 
4556   trc2d(ji,jj,91) = trc2d(ji,jj,91) + f_wind !! surface wind 
4557   trc2d(ji,jj,92) = trc2d(ji,jj,92) + f_pco2atm !! atmospheric pCO2 
4558   trc2d(ji,jj,93) = trc2d(ji,jj,93) + f_ph !! ocean pH 
4559   trc2d(ji,jj,94) = trc2d(ji,jj,94) + f_pco2w !! ocean pCO2 
4560   trc2d(ji,jj,95) = trc2d(ji,jj,95) + f_h2co3 !! ocean H2CO3 conc. 
4561   trc2d(ji,jj,96) = trc2d(ji,jj,96) + f_hco3 !! ocean HCO3 conc. 
4562   trc2d(ji,jj,97) = trc2d(ji,jj,97) + f_co3 !! ocean CO3 conc. 
4563   trc2d(ji,jj,98) = trc2d(ji,jj,98) + f_co2flux !! airsea CO2 flux 
4564   trc2d(ji,jj,99) = trc2d(ji,jj,99) + f_omcal(ji,jj) !! ocean omega calcite 
4565   trc2d(ji,jj,100) = trc2d(ji,jj,100) + f_omarg(ji,jj) !! ocean omega aragonite 
4566   trc2d(ji,jj,101) = trc2d(ji,jj,101) + f_TDIC !! ocean TDIC 
4567   trc2d(ji,jj,102) = trc2d(ji,jj,102) + f_TALK !! ocean TALK 
4568   trc2d(ji,jj,103) = trc2d(ji,jj,103) + f_kw660 !! surface kw660 
4569   trc2d(ji,jj,104) = trc2d(ji,jj,104) + f_pp0 !! surface pressure 
4570   trc2d(ji,jj,105) = trc2d(ji,jj,105) + f_o2flux !! airsea O2 flux 
4571   trc2d(ji,jj,106) = trc2d(ji,jj,106) + f_o2sat !! ocean O2 saturation 
4572   trc2d(ji,jj,107) = f2_ccd_cal(ji,jj) !! depth calcite CCD 
4573   trc2d(ji,jj,108) = f2_ccd_arg(ji,jj) !! depth aragonite CCD 
4574   endif 
4575   if (jk .eq. jmbathy) then 
4576   trc2d(ji,jj,109) = f3_omcal(ji,jj,jk) !! seafloor omega calcite 
4577   trc2d(ji,jj,110) = f3_omarg(ji,jj,jk) !! seafloor omega aragonite 
4578   endif 
4579   !! diagnostic fields 111 to 117 calculated below 
4580   if (jk.eq.i0100) trc2d(ji,jj,118) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 100 m 
4581   if (jk.eq.i0500) trc2d(ji,jj,119) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 500 m 
4582   if (jk.eq.i1000) trc2d(ji,jj,120) = ffastca(ji,jj)/MAX(ffastc(ji,jj), rsmall) !! rain ratio at 1000 m 
4583   !! AXY (18/01/12): benthic flux diagnostics 
4584   if (jk.eq.jmbathy) then 
4585   trc2d(ji,jj,121) = f_sbenin_n(ji,jj) + f_fbenin_n(ji,jj) 
4586   trc2d(ji,jj,122) = f_sbenin_fe(ji,jj) + f_fbenin_fe(ji,jj) 
4587   trc2d(ji,jj,123) = f_sbenin_c(ji,jj) + f_fbenin_c(ji,jj) 
4588   trc2d(ji,jj,124) = f_fbenin_si(ji,jj) 
4589   trc2d(ji,jj,125) = f_fbenin_ca(ji,jj) 
4590   trc2d(ji,jj,126) = f_benout_n(ji,jj) 
4591   trc2d(ji,jj,127) = f_benout_fe(ji,jj) 
4592   trc2d(ji,jj,128) = f_benout_c(ji,jj) 
4593   trc2d(ji,jj,129) = f_benout_si(ji,jj) 
4594   trc2d(ji,jj,130) = f_benout_ca(ji,jj) 
4595   endif 
4596   !! diagnostics fields 131 to 135 calculated below 
4597   trc2d(ji,jj,136) = f_runoff(ji,jj) 
4598   !! AXY (19/07/12): amended to allow for riverine nutrient addition below surface 
4599   trc2d(ji,jj,137) = trc2d(ji,jj,137) + (f_riv_loc_n * fthk) 
4600   trc2d(ji,jj,138) = trc2d(ji,jj,138) + (f_riv_loc_si * fthk) 
4601   trc2d(ji,jj,139) = trc2d(ji,jj,139) + (f_riv_loc_c * fthk) 
4602   trc2d(ji,jj,140) = trc2d(ji,jj,140) + (f_riv_loc_alk * fthk) 
4603   trc2d(ji,jj,141) = trc2d(ji,jj,141) + (fslowc * fthk) !! slow sinking detritus C production 
4604   if (jk.eq.i0100) trc2d(ji,jj,142) = fslowcflux(ji,jj) !! slow detritus flux at 100 m 
4605   if (jk.eq.i0200) trc2d(ji,jj,143) = fslowcflux(ji,jj) !! slow detritus flux at 200 m 
4606   if (jk.eq.i0500) trc2d(ji,jj,144) = fslowcflux(ji,jj) !! slow detritus flux at 500 m 
4607   if (jk.eq.i1000) trc2d(ji,jj,145) = fslowcflux(ji,jj) !! slow detritus flux at 1000 m 
4608   trc2d(ji,jj,146) = trc2d(ji,jj,146) + ftot_c(ji,jj) !! carbon inventory 
4609   trc2d(ji,jj,147) = trc2d(ji,jj,147) + ftot_a(ji,jj) !! alkalinity inventory 
4610   trc2d(ji,jj,148) = trc2d(ji,jj,148) + ftot_o2(ji,jj) !! oxygen inventory 
4611   if (jk.eq.jmbathy) then 
4612   trc2d(ji,jj,149) = f_benout_lyso_ca(ji,jj) 
4613   endif 
4614   trc2d(ji,jj,150) = fcomm_resp(ji,jj) * fthk !! community respiration 
4615   !! 
4616   !! AXY (14/02/14): a Valentines Day gift to BASIN  a shedload of new 
4617   !! diagnostics that they'll most likely never need! 
4618   !! (actually, as with all such gifts, I'm giving them 
4619   !! some things I'd like myself!) 
4620   !! 
4621   !!  
4622   !! linear losses 
4623   !! nondiatom 
4624   trc2d(ji,jj,151) = trc2d(ji,jj,151) + (fdpn2 * fthk) 
4625   !! diatom 
4626   trc2d(ji,jj,152) = trc2d(ji,jj,152) + (fdpd2 * fthk) 
4627   !! microzooplankton 
4628   trc2d(ji,jj,153) = trc2d(ji,jj,153) + (fdzmi2 * fthk) 
4629   !! mesozooplankton 
4630   trc2d(ji,jj,154) = trc2d(ji,jj,154) + (fdzme2 * fthk) 
4631   !!  
4632   !! microzooplankton grazing 
4633   !! microzooplankton messy > N 
4634   trc2d(ji,jj,155) = trc2d(ji,jj,155) + (xphi * (fgmipn + fgmid) * fthk) 
4635   !! microzooplankton messy > D 
4636   trc2d(ji,jj,156) = trc2d(ji,jj,156) + ((1.  xbetan) * finmi * fthk) 
4637   !! microzooplankton messy > DIC 
4638   trc2d(ji,jj,157) = trc2d(ji,jj,157) + (xphi * ((xthetapn * fgmipn) + fgmidc) * fthk) 
4639   !! microzooplankton messy > Dc 
4640   trc2d(ji,jj,158) = trc2d(ji,jj,158) + ((1.  xbetac) * ficmi * fthk) 
4641   !! microzooplankton excretion 
4642   trc2d(ji,jj,159) = trc2d(ji,jj,159) + (fmiexcr * fthk) 
4643   !! microzooplankton respiration 
4644   trc2d(ji,jj,160) = trc2d(ji,jj,160) + (fmiresp * fthk) 
4645   !! microzooplankton growth 
4646   trc2d(ji,jj,161) = trc2d(ji,jj,161) + (fmigrow * fthk) 
4647   !!  
4648   !! mesozooplankton grazing 
4649   !! mesozooplankton messy > N 
4650   trc2d(ji,jj,162) = trc2d(ji,jj,162) + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed) * fthk) 
4651   !! mesozooplankton messy > D 
4652   trc2d(ji,jj,163) = trc2d(ji,jj,163) + ((1.  xbetan) * finme * fthk) 
4653   !! mesozooplankton messy > DIC 
4654   trc2d(ji,jj,164) = trc2d(ji,jj,164) + (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & 
4655   & (xthetazmi * fgmezmi) + fgmedc) * fthk) 
4656   !! mesozooplankton messy > Dc 
4657   trc2d(ji,jj,165) = trc2d(ji,jj,165) + ((1.  xbetac) * ficme * fthk) 
4658   !! mesozooplankton excretion 
4659   trc2d(ji,jj,166) = trc2d(ji,jj,166) + (fmeexcr * fthk) 
4660   !! mesozooplankton respiration 
4661   trc2d(ji,jj,167) = trc2d(ji,jj,167) + (fmeresp * fthk) 
4662   !! mesozooplankton growth 
4663   trc2d(ji,jj,168) = trc2d(ji,jj,168) + (fmegrow * fthk) 
4664   !!  
4665   !! miscellaneous 
4666   trc2d(ji,jj,169) = trc2d(ji,jj,169) + (fddc * fthk) !! detrital C remineralisation 
4667   trc2d(ji,jj,170) = trc2d(ji,jj,170) + (fgmidc * fthk) !! microzoo grazing on detrital carbon 
4668   trc2d(ji,jj,171) = trc2d(ji,jj,171) + (fgmedc * fthk) !! mesozoo grazing on detrital carbon 
4669   !! 
4670   !!  
4671   !! 
4672   !! AXY (23/10/14): extract primary production related surface fields to 
4673   !! deal with diel cycle issues; hijacking BASIN 150m 
4674   !! diagnostics to do so (see commented out diagnostics 
4675   !! below this section) 
4676   !! 
4677   !! extract fields at surface 
4678   !! if (jk .eq. 1) then 
4679   !! trc2d(ji,jj,172) = zchn !! Pn chlorophyll 
4680   !! trc2d(ji,jj,173) = zphn !! Pn biomass 
4681   !! trc2d(ji,jj,174) = fjln !! Pn Jterm 
4682   !! trc2d(ji,jj,175) = (fprn * zphn) !! Pn PP 
4683   !! trc2d(ji,jj,176) = zchd !! Pd chlorophyll 
4684   !! trc2d(ji,jj,177) = zphd !! Pd biomass 
4685   !! trc2d(ji,jj,178) = fjld !! Pd Jterm 
4686   !! trc2d(ji,jj,179) = xpar(ji,jj,jk) !! Pd PP 
4687   !! trc2d(ji,jj,180) = loc_T !! local temperature 
4688   !! endif 
4689   !! !! 
4690   !! !! extract fields at 50m (actually 4450m) 
4691   !! if (jk .eq. 18) then 
4692   !! trc2d(ji,jj,181) = zchn !! Pn chlorophyll 
4693   !! trc2d(ji,jj,182) = zphn !! Pn biomass 
4694   !! trc2d(ji,jj,183) = fjln !! Pn Jterm 
4695   !! trc2d(ji,jj,184) = (fprn * zphn) !! Pn PP 
4696   !! trc2d(ji,jj,185) = zchd !! Pd chlorophyll 
4697   !! trc2d(ji,jj,186) = zphd !! Pd biomass 
4698   !! trc2d(ji,jj,187) = fjld !! Pd Jterm 
4699   !! trc2d(ji,jj,188) = xpar(ji,jj,jk) !! Pd PP 
4700   !! trc2d(ji,jj,189) = loc_T !! local temperature 
4701   !! endif 
4702   !! !! 
4703   !! !! extract fields at 100m 
4704   !! if (jk .eq. i0100) then 
4705   !! trc2d(ji,jj,190) = zchn !! Pn chlorophyll 
4706   !! trc2d(ji,jj,191) = zphn !! Pn biomass 
4707   !! trc2d(ji,jj,192) = fjln !! Pn Jterm 
4708   !! trc2d(ji,jj,193) = (fprn * zphn) !! Pn PP 
4709   !! trc2d(ji,jj,194) = zchd !! Pd chlorophyll 
4710   !! trc2d(ji,jj,195) = zphd !! Pd biomass 
4711   !! trc2d(ji,jj,196) = fjld !! Pd Jterm 
4712   !! trc2d(ji,jj,197) = xpar(ji,jj,jk) !! Pd PP 
4713   !! trc2d(ji,jj,198) = loc_T !! local temperature 
4714   !! endif 
4715   !! 
4716   !! extract relevant BASIN fields at 150m 
4717   if (jk .eq. i0150) then 
4718   trc2d(ji,jj,172) = trc2d(ji,jj,4) !! Pn PP 
4719   trc2d(ji,jj,173) = trc2d(ji,jj,151) !! Pn linear loss 
4720   trc2d(ji,jj,174) = trc2d(ji,jj,5) !! Pn nonlinear loss 
