403   !! 
404   
405   # if defined key_debug_medusa 
406   IF (lwp) write (numout,*) 'trc_bio_medusa: variables defined' 
407   CALL flush(numout) 
408   # endif 
409   
410   !! AXY (20/11/14): alter this to report on first MEDUSA call 
411   !! IF( kt == nit000 ) THEN 
412   IF( kt == nittrc000 ) THEN 
413   IF(lwp) WRITE(numout,*) 
414   IF(lwp) WRITE(numout,*) ' trc_bio: MEDUSA biomodel' 
415   IF(lwp) WRITE(numout,*) ' ~~~~~~~' 
416   IF(lwp) WRITE(numout,*) ' kt =',kt 
417   ENDIF 
418   
419   !! AXY (13/01/12): is benthic model properly interactive? 0 = no, 1 = yes 
420   ibenthic = 1 
421   
422   !! not sure what this is for; it's not used anywhere; commenting out 
423   !! fbodn(:,:) = 0.e0 
424   
425   !! 
426   IF( ln_diatrc ) THEN 
427   !! blank 2D diagnostic array 
428   trc2d(:,:,:) = 0.e0 
429   !! 
430   !! blank 3D diagnostic array 
431   trc3d(:,:,:,:) = 0.e0 
432   ENDIF 
433   
434   !! 
435   !! b0 is present for debugging purposes; using b0 = 0 sets the tendency 
436   !! terms of all biological equations to 0. 
437   !! 
438   !! 
439   !! AXY (03/09/14): probably not the smartest move ever, but it'll fit 
440   !! the bill for now; another item on the thingstosort 
441   !! outinthefuture list ... 
442   # if defined key_kill_medusa 
443   b0 = 0. 
444   # else 
445   b0 = 1. 
446   # endif 
447   !! 
448   !! fast detritus ballast scheme (0 = no; 1 = yes) 
449   !! alternative to ballast scheme is same scheme but with no ballast 
450   !! protection (not dissimilar to Martin et al., 1987) 
451   !! 
452   !! 
453   iball = 1 
454   
455   !! 
456   !! full flux diagnostics (0 = no; 1 = yes); appear in ocean.output 
457   !! these should *only* be used in 1D since they give comprehensive 
458   !! output for ecological functions in the model; primarily used in 
459   !! debugging 
460   !! 
461   !! 
462   idf = 0 
463   !! 
464   !! timer mechanism 
465   if (kt/120*120.eq.kt) then 
466   idfval = 1 
467   else 
468   idfval = 0 
469   endif 
470   
471   !! 
472   !! blank fastsinking detritus 2D fields 
473   !! 
474   !! 
475   ffastn(:,:) = 0.0 !! organic nitrogen 
476   ffastsi(:,:) = 0.0 !! biogenic silicon 
477   ffastfe(:,:) = 0.0 !! organic iron 
478   ffastc(:,:) = 0.0 !! organic carbon 
479   ffastca(:,:) = 0.0 !! biogenic calcium carbonate 
480   !! 
481   fsedn(:,:) = 0.0 !! Seafloor flux of N 
482   fsedsi(:,:) = 0.0 !! Seafloor flux of Si 
483   fsedfe(:,:) = 0.0 !! Seafloor flux of Fe 
484   fsedc(:,:) = 0.0 !! Seafloor flux of C 
485   fsedca(:,:) = 0.0 !! Seafloor flux of CaCO3 
486   !! 
487   fregenfast(:,:) = 0.0 !! integrated N regeneration (fast detritus) 
488   fregenfastsi(:,:) = 0.0 !! integrated Si regeneration (fast detritus) 
 516  ftot_n(:,:) = 0.0 !! N inventory 
 517  ftot_si(:,:) = 0.0 !! Si inventory 
 518  ftot_fe(:,:) = 0.0 !! Fe inventory 
 519  ftot_pn(:,:) = 0.0 !! integrated nondiatom phytoplankton 
 520  ftot_pd(:,:) = 0.0 !! integrated diatom phytoplankton 
 521  ftot_zmi(:,:) = 0.0 !! integrated microzooplankton 
 522  ftot_zme(:,:) = 0.0 !! integrated mesozooplankton 
 523  ftot_det(:,:) = 0.0 !! integrated slow detritus, nitrogen 
 524  ftot_dtc(:,:) = 0.0 !! integrated slow detritus, carbon 
 525  !! 
 526  fzmi_i(:,:) = 0.0 !! material grazed by microzooplankton 
 527  fzmi_o(:,:) = 0.0 !! ... sum of fate of this material 
 528  fzme_i(:,:) = 0.0 !! material grazed by mesozooplankton 
 529  fzme_o(:,:) = 0.0 !! ... sum of fate of this material 
 530  !! 
 531  f_sbenin_n(:,:) = 0.0 !! slow detritus N > benthic pool 
 532  f_sbenin_fe(:,:) = 0.0 !! slow detritus Fe > benthic pool 
 533  f_sbenin_c(:,:) = 0.0 !! slow detritus C > benthic pool 
 534  f_fbenin_n(:,:) = 0.0 !! fast detritus N > benthic pool 
 535  f_fbenin_fe(:,:) = 0.0 !! fast detritus Fe > benthic pool 
 536  f_fbenin_si(:,:) = 0.0 !! fast detritus Si > benthic pool 
 537  f_fbenin_c(:,:) = 0.0 !! fast detritus C > benthic pool 
 538  f_fbenin_ca(:,:) = 0.0 !! fast detritus Ca > benthic pool 
 539  !! 
 540  f_benout_n(:,:) = 0.0 !! benthic N pool > dissolved 
 541  f_benout_fe(:,:) = 0.0 !! benthic Fe pool > dissolved 
 542  f_benout_si(:,:) = 0.0 !! benthic Si pool > dissolved 
 543  f_benout_c(:,:) = 0.0 !! benthic C pool > dissolved 
 544  f_benout_ca(:,:) = 0.0 !! benthic Ca pool > dissolved 
 545  !! 
 546  f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool > dissolved (when it shouldn't!) 
 547  !! 
 548  f_runoff(:,:) = 0.0 !! riverine runoff 
 549  f_riv_n(:,:) = 0.0 !! riverine N input 
 550  f_riv_si(:,:) = 0.0 !! riverine Si input 
 551  f_riv_c(:,:) = 0.0 !! riverine C input 
 552  f_riv_alk(:,:) = 0.0 !! riverine alk input 
 553  !! 
 554  !! Jpalm  06032017  Forgotten var to init 
 555  f_omarg(:,:) = 0.0 !! 
 556  f_omcal(:,:) = 0.0 
 557  xFree(:,:) = 0.0 !! state variables for ironligand system 
 558  fcomm_resp(:,:) = 0.0 
 559  fprn_ml(:,:) = 0.0 !! mixed layer PP diagnostics 
 560  fprd_ml(:,:) = 0.0 !! mixed layer PP diagnostics 
 561  
 562  !! 
 563  !! allocate and initiate 2D diag 
 564  !! 
 565  !! 
 566  IF ( lk_iomput .AND. .NOT. ln_diatrc ) THEN 
 567  !! Juju :: add kt condition !! 
 568  if ( kt == nittrc000 ) CALL trc_nam_iom_medusa !! initialise iom_use test 
 569  !! 
 570  CALL wrk_alloc( jpi, jpj, zw2d ) 
 571  zw2d(:,:) = 0.0 !! 
 572  IF ( med_diag%PRN%dgsave ) THEN 
 573  CALL wrk_alloc( jpi, jpj, fprn2d ) 
 574  fprn2d(:,:) = 0.0 !! 
 575  ENDIF 
 576  IF ( med_diag%MPN%dgsave ) THEN 
 577  CALL wrk_alloc( jpi, jpj, fdpn2d ) 
 578  fdpn2d(:,:) = 0.0 !! 
 579  ENDIF 
 580  IF ( med_diag%PRD%dgsave ) THEN 
 581  CALL wrk_alloc( jpi, jpj, fprd2d ) 
 582  fprd2d(:,:) = 0.0 !! 
 583  ENDIF 
 584  IF( med_diag%MPD%dgsave ) THEN 
 585  CALL wrk_alloc( jpi, jpj, fdpd2d ) 
 586  fdpd2d(:,:) = 0.0 !! 
 587  ENDIF 
 588  IF( med_diag%OPAL%dgsave ) THEN 
 589  CALL wrk_alloc( jpi, jpj, fprds2d ) 
 590  fprds2d(:,:) = 0.0 !! 
 591  ENDIF 
 592  IF( med_diag%OPALDISS%dgsave ) THEN 
 593  CALL wrk_alloc( jpi, jpj, fsdiss2d ) 
 594  fsdiss2d(:,:) = 0.0 !! 
 595  ENDIF 
 596  IF( med_diag%GMIPn%dgsave ) THEN 
 597  CALL wrk_alloc( jpi, jpj, fgmipn2d ) 
 598  fgmipn2d(:,:) = 0.0 !! 
 599  ENDIF 
 600  IF( med_diag%GMID%dgsave ) THEN 
 601  CALL wrk_alloc( jpi, jpj, fgmid2d ) 
 602  fgmid2d(:,:) = 0.0 !! 
 603  ENDIF 
 604  IF( med_diag%MZMI%dgsave ) THEN 
 605  CALL wrk_alloc( jpi, jpj, fdzmi2d ) 
 606  fdzmi2d(:,:) = 0.0 !! 
 607  ENDIF 
 608  IF( med_diag%GMEPN%dgsave ) THEN 
 609  CALL wrk_alloc( jpi, jpj, fgmepn2d ) 
 610  fgmepn2d(:,:) = 0.0 !! 
 611  ENDIF 
 612  IF( med_diag%GMEPD%dgsave ) THEN 
 613  CALL wrk_alloc( jpi, jpj, fgmepd2d ) 
 614  fgmepd2d(:,:) = 0.0 !! 
 615  ENDIF 
 616  IF( med_diag%GMEZMI%dgsave ) THEN 
 617  CALL wrk_alloc( jpi, jpj, fgmezmi2d ) 
 618  fgmezmi2d(:,:) = 0.0 !! 
 619  ENDIF 
 620  IF( med_diag%GMED%dgsave ) THEN 
 621  CALL wrk_alloc( jpi, jpj, fgmed2d ) 
 622  fgmed2d(:,:) = 0.0 !! 
 623  ENDIF 
 624  IF( med_diag%MZME%dgsave ) THEN 
 625  CALL wrk_alloc( jpi, jpj, fdzme2d ) 
 626  fdzme2d(:,:) = 0.0 !! 
 627  ENDIF 
 628  IF( med_diag%DETN%dgsave ) THEN 
 629  CALL wrk_alloc( jpi, jpj, fslown2d ) 
 630  fslown2d(:,:) = 0.0 !! 
 631  ENDIF 
 632  IF( med_diag%MDET%dgsave ) THEN 
 633  CALL wrk_alloc( jpi, jpj, fdd2d ) 
 634  fdd2d(:,:) = 0.0 !! 
 635  ENDIF 
 636  IF( med_diag%AEOLIAN%dgsave ) THEN 
 637  CALL wrk_alloc( jpi, jpj, ffetop2d ) 
 638  ffetop2d(:,:) = 0.0 !! 
 639  ENDIF 
 640  IF( med_diag%BENTHIC%dgsave ) THEN 
 641  CALL wrk_alloc( jpi, jpj, ffebot2d ) 
 642  ffebot2d(:,:) = 0.0 !! 
 643  ENDIF 
 644  IF( med_diag%SCAVENGE%dgsave ) THEN 
 645  CALL wrk_alloc( jpi, jpj, ffescav2d ) 
 646  ffescav2d(:,:) = 0.0 !! 
 647  ENDIF 
 648  IF( med_diag%PN_JLIM%dgsave ) THEN 
 649  CALL wrk_alloc( jpi, jpj, fjln2d ) 
 650  fjln2d(:,:) = 0.0 !! 
 651  ENDIF 
 652  IF( med_diag%PN_NLIM%dgsave ) THEN 
 653  CALL wrk_alloc( jpi, jpj, fnln2d ) 
 654  fnln2d(:,:) = 0.0 !! 
 655  ENDIF 
 656  IF( med_diag%PN_FELIM%dgsave ) THEN 
 657  CALL wrk_alloc( jpi, jpj, ffln2d ) 
 658  ffln2d(:,:) = 0.0 !! 
 659  ENDIF 
 660  IF( med_diag%PD_JLIM%dgsave ) THEN 
 661  CALL wrk_alloc( jpi, jpj, fjld2d ) 
 662  fjld2d(:,:) = 0.0 !! 
 663  ENDIF 
 664  IF( med_diag%PD_NLIM%dgsave ) THEN 
 665  CALL wrk_alloc( jpi, jpj, fnld2d ) 
 666  fnld2d(:,:) = 0.0 !! 
 667  ENDIF 
 668  IF( med_diag%PD_FELIM%dgsave ) THEN 
 669  CALL wrk_alloc( jpi, jpj, ffld2d ) 
 670  ffld2d(:,:) = 0.0 !! 
 671  ENDIF 
 672  IF( med_diag%PD_SILIM%dgsave ) THEN 
 673  CALL wrk_alloc( jpi, jpj, fsld2d2 ) 
 674  fsld2d2(:,:) = 0.0 !! 
 675  ENDIF 
 676  IF( med_diag%PDSILIM2%dgsave ) THEN 
 677  CALL wrk_alloc( jpi, jpj, fsld2d ) 
 678  fsld2d(:,:) = 0.0 !! 
 679  ENDIF 
 680  !! 
 681  !! skip SDT_XXXX diagnostics here 
 682  !! 
 683  IF( med_diag%TOTREG_N%dgsave ) THEN 
 684  CALL wrk_alloc( jpi, jpj, fregen2d ) 
 685  fregen2d(:,:) = 0.0 !! 
 686  ENDIF 
 687  IF( med_diag%TOTRG_SI%dgsave ) THEN 
 688  CALL wrk_alloc( jpi, jpj, fregensi2d ) 
 689  fregensi2d(:,:) = 0.0 !! 
 690  ENDIF 
 691  !! 
 692  !! skip REG_XXXX diagnostics here 
 693  !! 
 694  IF( med_diag%FASTN%dgsave ) THEN 
 695  CALL wrk_alloc( jpi, jpj, ftempn2d ) 
 696  ftempn2d(:,:) = 0.0 !! 
 697  ENDIF 
 698  IF( med_diag%FASTSI%dgsave ) THEN 
 699  CALL wrk_alloc( jpi, jpj, ftempsi2d ) 
 700  ftempsi2d(:,:) = 0.0 !! 
 701  ENDIF 
 702  IF( med_diag%FASTFE%dgsave ) THEN 
 703  CALL wrk_alloc( jpi, jpj, ftempfe2d ) 
 704  ftempfe2d(:,:) = 0.0 !! 
 705  ENDIF 
 706  IF( med_diag%FASTC%dgsave ) THEN 
 707  CALL wrk_alloc( jpi, jpj, ftempc2d ) 
 708  ftempc2d(:,:) = 0.0 !! 
 709  ENDIF 
 710  IF( med_diag%FASTCA%dgsave ) THEN 
 711  CALL wrk_alloc( jpi, jpj, ftempca2d ) 
 712  ftempca2d(:,:) = 0.0 !! 
 713  ENDIF 
 714  !! 
 715  !! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here 
 716  !! 
 717  IF( med_diag%REMINN%dgsave ) THEN 
 718  CALL wrk_alloc( jpi, jpj, freminn2d ) 
 719  freminn2d(:,:) = 0.0 !! 
 720  ENDIF 
 721  IF( med_diag%REMINSI%dgsave ) THEN 
 722  CALL wrk_alloc( jpi, jpj, freminsi2d ) 
 723  freminsi2d(:,:) = 0.0 !! 
 724  ENDIF 
 725  IF( med_diag%REMINFE%dgsave ) THEN 
 726  CALL wrk_alloc( jpi, jpj, freminfe2d ) 
 727  freminfe2d(:,:) = 0.0 !! 
 728  ENDIF 
 729  IF( med_diag%REMINC%dgsave ) THEN 
 730  CALL wrk_alloc( jpi, jpj, freminc2d ) 
 731  freminc2d(:,:) = 0.0 !! 
 732  ENDIF 
 733  IF( med_diag%REMINCA%dgsave ) THEN 
 734  CALL wrk_alloc( jpi, jpj, freminca2d ) 
 735  freminca2d(:,:) = 0.0 !! 
 736  ENDIF 
490   fregenfastc(:,:) = 0.0 !! integrated C regeneration (fast detritus) 
491   # endif 
492   !! 
493   fccd(:,:) = 0.0 !! last depth level before CCD 
494   
495   !! 
496   !! blank nutrient/flux inventories 
497   !! 
498   !! 
499   fflx_n(:,:) = 0.0 !! nitrogen flux total 
500   fflx_si(:,:) = 0.0 !! silicon flux total 
501   fflx_fe(:,:) = 0.0 !! iron flux total 
502   fifd_n(:,:) = 0.0 !! nitrogen fast detritus production 
503   fifd_si(:,:) = 0.0 !! silicon fast detritus production 
504   fifd_fe(:,:) = 0.0 !! iron fast detritus production 
505   fofd_n(:,:) = 0.0 !! nitrogen fast detritus remineralisation 
506   fofd_si(:,:) = 0.0 !! silicon fast detritus remineralisation 
507   fofd_fe(:,:) = 0.0 !! iron fast detritus remineralisation 
508   # if defined key_roam 
509   fflx_c(:,:) = 0.0 !! carbon flux total 
510   fflx_a(:,:) = 0.0 !! alkalinity flux total 
511   fflx_o2(:,:) = 0.0 !! oxygen flux total 
512   ftot_c(:,:) = 0.0 !! carbon inventory 
513   ftot_a(:,:) = 0.0 !! alkalinity inventory 
514   ftot_o2(:,:) = 0.0 !! oxygen inventory 
515   fifd_c(:,:) = 0.0 !! carbon fast detritus production 
516   fifd_a(:,:) = 0.0 !! alkalinity fast detritus production 
517   fifd_o2(:,:) = 0.0 !! oxygen fast detritus production 
518   fofd_c(:,:) = 0.0 !! carbon fast detritus remineralisation 
519   fofd_a(:,:) = 0.0 !! alkalinity fast detritus remineralisation 
520   fofd_o2(:,:) = 0.0 !! oxygen fast detritus remineralisation 
521   !! 
522   fnit_prod(:,:) = 0.0 !! (organic) nitrogen production 
523   fnit_cons(:,:) = 0.0 !! (organic) nitrogen consumption 
524   fsil_prod(:,:) = 0.0 !! (inorganic) silicon production 
525   fsil_cons(:,:) = 0.0 !! (inorganic) silicon consumption 
526   fcar_prod(:,:) = 0.0 !! (organic) carbon production 
527   fcar_cons(:,:) = 0.0 !! (organic) carbon consumption 
528   !! 
529   foxy_prod(:,:) = 0.0 !! oxygen production 
530   foxy_cons(:,:) = 0.0 !! oxygen consumption 
531   foxy_anox(:,:) = 0.0 !! unrealised oxygen consumption 
532   !! 
533   # endif 
534   ftot_n(:,:) = 0.0 !! N inventory 
535   ftot_si(:,:) = 0.0 !! Si inventory 
536   ftot_fe(:,:) = 0.0 !! Fe inventory 
537   ftot_pn(:,:) = 0.0 !! integrated nondiatom phytoplankton 
538   ftot_pd(:,:) = 0.0 !! integrated diatom phytoplankton 
539   ftot_zmi(:,:) = 0.0 !! integrated microzooplankton 
540   ftot_zme(:,:) = 0.0 !! integrated mesozooplankton 
541   ftot_det(:,:) = 0.0 !! integrated slow detritus, nitrogen 
542   ftot_dtc(:,:) = 0.0 !! integrated slow detritus, carbon 
543   !! 
544   fzmi_i(:,:) = 0.0 !! material grazed by microzooplankton 
545   fzmi_o(:,:) = 0.0 !! ... sum of fate of this material 
546   fzme_i(:,:) = 0.0 !! material grazed by mesozooplankton 
547   fzme_o(:,:) = 0.0 !! ... sum of fate of this material 
548   !! 
549   f_sbenin_n(:,:) = 0.0 !! slow detritus N > benthic pool 
550   f_sbenin_fe(:,:) = 0.0 !! slow detritus Fe > benthic pool 
551   f_sbenin_c(:,:) = 0.0 !! slow detritus C > benthic pool 
552   f_fbenin_n(:,:) = 0.0 !! fast detritus N > benthic pool 
553   f_fbenin_fe(:,:) = 0.0 !! fast detritus Fe > benthic pool 
554   f_fbenin_si(:,:) = 0.0 !! fast detritus Si > benthic pool 
555   f_fbenin_c(:,:) = 0.0 !! fast detritus C > benthic pool 
556   f_fbenin_ca(:,:) = 0.0 !! fast detritus Ca > benthic pool 
557   !! 
558   f_benout_n(:,:) = 0.0 !! benthic N pool > dissolved 
559   f_benout_fe(:,:) = 0.0 !! benthic Fe pool > dissolved 
560   f_benout_si(:,:) = 0.0 !! benthic Si pool > dissolved 
561   f_benout_c(:,:) = 0.0 !! benthic C pool > dissolved 
562   f_benout_ca(:,:) = 0.0 !! benthic Ca pool > dissolved 
563   !! 
564   f_benout_lyso_ca(:,:) = 0.0 !! benthic Ca pool > dissolved (when it shouldn't!) 
565   !! 
566   f_runoff(:,:) = 0.0 !! riverine runoff 
567   f_riv_n(:,:) = 0.0 !! riverine N input 
568   f_riv_si(:,:) = 0.0 !! riverine Si input 
569   f_riv_c(:,:) = 0.0 !! riverine C input 
570   f_riv_alk(:,:) = 0.0 !! riverine alk input 
571   !! 
572   !! Jpalm  06032017  Forgotten var to init 
573   f_omarg(:,:) = 0.0 !! 
574   f_omcal(:,:) = 0.0 
575   xFree(:,:) = 0.0 !! state variables for ironligand system 
576   fcomm_resp(:,:) = 0.0 
577   fprn_ml(:,:) = 0.0 !! mixed layer PP diagnostics 
578   fprd_ml(:,:) = 0.0 !! mixed layer PP diagnostics 
579   
580   !! 
581   !! allocate and initiate 2D diag 
582   !!  
583   !! Juju :: add kt condition !! 
584   IF ( lk_iomput .AND. .NOT. ln_diatrc ) THEN 
585   !! 
586   if ( kt == nittrc000 ) CALL trc_nam_iom_medusa !! initialise iom_use test 
587   !! 
588   CALL wrk_alloc( jpi, jpj, zw2d ) 
589   zw2d(:,:) = 0.0 !! 
590   IF ( med_diag%PRN%dgsave ) THEN 
591   CALL wrk_alloc( jpi, jpj, fprn2d ) 
592   fprn2d(:,:) = 0.0 !! 
593   ENDIF 
594   IF ( med_diag%MPN%dgsave ) THEN 
595   CALL wrk_alloc( jpi, jpj, fdpn2d ) 
596   fdpn2d(:,:) = 0.0 !! 
597   ENDIF 
598   IF ( med_diag%PRD%dgsave ) THEN 
599   CALL wrk_alloc( jpi, jpj, fprd2d ) 
600   fprd2d(:,:) = 0.0 !! 
601   ENDIF 
602   IF( med_diag%MPD%dgsave ) THEN 
603   CALL wrk_alloc( jpi, jpj, fdpd2d ) 
604   fdpd2d(:,:) = 0.0 !! 
605   ENDIF 
606   IF( med_diag%OPAL%dgsave ) THEN 
607   CALL wrk_alloc( jpi, jpj, fprds2d ) 
608   fprds2d(:,:) = 0.0 !! 
609   ENDIF 
610   IF( med_diag%OPALDISS%dgsave ) THEN 
611   CALL wrk_alloc( jpi, jpj, fsdiss2d ) 
612   fsdiss2d(:,:) = 0.0 !! 
613   ENDIF 
614   IF( med_diag%GMIPn%dgsave ) THEN 
615   CALL wrk_alloc( jpi, jpj, fgmipn2d ) 
616   fgmipn2d(:,:) = 0.0 !! 
617   ENDIF 
618   IF( med_diag%GMID%dgsave ) THEN 
619   CALL wrk_alloc( jpi, jpj, fgmid2d ) 
620   fgmid2d(:,:) = 0.0 !! 
621   ENDIF 
622   IF( med_diag%MZMI%dgsave ) THEN 
623   CALL wrk_alloc( jpi, jpj, fdzmi2d ) 
624   fdzmi2d(:,:) = 0.0 !! 
625   ENDIF 
626   IF( med_diag%GMEPN%dgsave ) THEN 
627   CALL wrk_alloc( jpi, jpj, fgmepn2d ) 
628   fgmepn2d(:,:) = 0.0 !! 
629   ENDIF 
630   IF( med_diag%GMEPD%dgsave ) THEN 
631   CALL wrk_alloc( jpi, jpj, fgmepd2d ) 
632   fgmepd2d(:,:) = 0.0 !! 
633   ENDIF 
634   IF( med_diag%GMEZMI%dgsave ) THEN 
635   CALL wrk_alloc( jpi, jpj, fgmezmi2d ) 
636   fgmezmi2d(:,:) = 0.0 !! 
637   ENDIF 
638   IF( med_diag%GMED%dgsave ) THEN 
639   CALL wrk_alloc( jpi, jpj, fgmed2d ) 
640   fgmed2d(:,:) = 0.0 !! 
641   ENDIF 
642   IF( med_diag%MZME%dgsave ) THEN 
643   CALL wrk_alloc( jpi, jpj, fdzme2d ) 
644   fdzme2d(:,:) = 0.0 !! 
645   ENDIF 
646   IF( med_diag%DETN%dgsave ) THEN 
647   CALL wrk_alloc( jpi, jpj, fslown2d ) 
648   fslown2d(:,:) = 0.0 !! 
649   ENDIF 
650   IF( med_diag%MDET%dgsave ) THEN 
651   CALL wrk_alloc( jpi, jpj, fdd2d ) 
652   fdd2d(:,:) = 0.0 !! 
653   ENDIF 
654   IF( med_diag%AEOLIAN%dgsave ) THEN 
655   CALL wrk_alloc( jpi, jpj, ffetop2d ) 
656   ffetop2d(:,:) = 0.0 !! 
657   ENDIF 
658   IF( med_diag%BENTHIC%dgsave ) THEN 
659   CALL wrk_alloc( jpi, jpj, ffebot2d ) 
660   ffebot2d(:,:) = 0.0 !! 
661   ENDIF 
662   IF( med_diag%SCAVENGE%dgsave ) THEN 
663   CALL wrk_alloc( jpi, jpj, ffescav2d ) 
664   ffescav2d(:,:) = 0.0 !! 
665   ENDIF 
666   IF( med_diag%PN_JLIM%dgsave ) THEN 
667   CALL wrk_alloc( jpi, jpj, fjln2d ) 
668   fjln2d(:,:) = 0.0 !! 
669   ENDIF 
670   IF( med_diag%PN_NLIM%dgsave ) THEN 
671   CALL wrk_alloc( jpi, jpj, fnln2d ) 
672   fnln2d(:,:) = 0.0 !! 
673   ENDIF 
674   IF( med_diag%PN_FELIM%dgsave ) THEN 
675   CALL wrk_alloc( jpi, jpj, ffln2d ) 
676   ffln2d(:,:) = 0.0 !! 
677   ENDIF 
678   IF( med_diag%PD_JLIM%dgsave ) THEN 
679   CALL wrk_alloc( jpi, jpj, fjld2d ) 
680   fjld2d(:,:) = 0.0 !! 
681   ENDIF 
682   IF( med_diag%PD_NLIM%dgsave ) THEN 
683   CALL wrk_alloc( jpi, jpj, fnld2d ) 
684   fnld2d(:,:) = 0.0 !! 
685   ENDIF 
686   IF( med_diag%PD_FELIM%dgsave ) THEN 
687   CALL wrk_alloc( jpi, jpj, ffld2d ) 
688   ffld2d(:,:) = 0.0 !! 
689   ENDIF 
690   IF( med_diag%PD_SILIM%dgsave ) THEN 
691   CALL wrk_alloc( jpi, jpj, fsld2d2 ) 
692   fsld2d2(:,:) = 0.0 !! 
693   ENDIF 
694   IF( med_diag%PDSILIM2%dgsave ) THEN 
695   CALL wrk_alloc( jpi, jpj, fsld2d ) 
696   fsld2d(:,:) = 0.0 !! 
697   ENDIF 
698   !! 
699   !! skip SDT_XXXX diagnostics here 
700   !! 
701   IF( med_diag%TOTREG_N%dgsave ) THEN 
702   CALL wrk_alloc( jpi, jpj, fregen2d ) 
703   fregen2d(:,:) = 0.0 !! 
704   ENDIF 
705   IF( med_diag%TOTRG_SI%dgsave ) THEN 
706   CALL wrk_alloc( jpi, jpj, fregensi2d ) 
707   fregensi2d(:,:) = 0.0 !! 
708   ENDIF 
709   !! 
710   !! skip REG_XXXX diagnostics here 
711   !! 
712   IF( med_diag%FASTN%dgsave ) THEN 
713   CALL wrk_alloc( jpi, jpj, ftempn2d ) 
714   ftempn2d(:,:) = 0.0 !! 
715   ENDIF 
716   IF( med_diag%FASTSI%dgsave ) THEN 
717   CALL wrk_alloc( jpi, jpj, ftempsi2d ) 
718   ftempsi2d(:,:) = 0.0 !! 
719   ENDIF 
720   IF( med_diag%FASTFE%dgsave ) THEN 
721   CALL wrk_alloc( jpi, jpj, ftempfe2d ) 
722   ftempfe2d(:,:) = 0.0 !! 
723   ENDIF 
724   IF( med_diag%FASTC%dgsave ) THEN 
725   CALL wrk_alloc( jpi, jpj, ftempc2d ) 
726   ftempc2d(:,:) = 0.0 !! 
727   ENDIF 
728   IF( med_diag%FASTCA%dgsave ) THEN 
729   CALL wrk_alloc( jpi, jpj, ftempca2d ) 
730   ftempca2d(:,:) = 0.0 !! 
731   ENDIF 
732   !! 
733   !! skip FDT_XXXX, RG_XXXXF, FDS_XXXX, RGS_XXXXF diagnostics here 
734   !! 
735   IF( med_diag%REMINN%dgsave ) THEN 
736   CALL wrk_alloc( jpi, jpj, freminn2d ) 
737   freminn2d(:,:) = 0.0 !! 
738   ENDIF 
739   IF( med_diag%REMINSI%dgsave ) THEN 
740   CALL wrk_alloc( jpi, jpj, freminsi2d ) 
741   freminsi2d(:,:) = 0.0 !! 
742   ENDIF 
743   IF( med_diag%REMINFE%dgsave ) THEN 
744   CALL wrk_alloc( jpi, jpj, freminfe2d ) 
745   freminfe2d(:,:) = 0.0 !! 
746   ENDIF 
747   IF( med_diag%REMINC%dgsave ) THEN 
748   CALL wrk_alloc( jpi, jpj, freminc2d ) 
749   freminc2d(:,:) = 0.0 !! 
750   ENDIF 
751   IF( med_diag%REMINCA%dgsave ) THEN 
752   CALL wrk_alloc( jpi, jpj, freminca2d ) 
753   freminca2d(:,:) = 0.0 !! 
754   ENDIF 
755   # if defined key_roam 
756   !! 
757   !! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here 
758   !! 
759   IF( med_diag%RR_0100%dgsave ) THEN 
760   CALL wrk_alloc( jpi, jpj, ffastca2d ) 
761   ffastca2d(:,:) = 0.0 !! 
762   ENDIF 
763   
764   IF( med_diag%ATM_PCO2%dgsave ) THEN 
765   CALL wrk_alloc( jpi, jpj, f_pco2a2d ) 
766   f_pco2a2d(:,:) = 0.0 !! 
767   ENDIF 
768   !! 
769   !! skip OCN_PH diagnostic here 
770   !! 
771   IF( med_diag%OCN_PCO2%dgsave ) THEN 
772   CALL wrk_alloc( jpi, jpj, f_pco2w2d ) 
773   f_pco2w2d(:,:) = 0.0 !! 
774   ENDIF 
775   !! 
776   !! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here 
777   !! 
778   IF( med_diag%CO2FLUX%dgsave ) THEN 
779   CALL wrk_alloc( jpi, jpj, f_co2flux2d ) 
780   f_co2flux2d(:,:) = 0.0 !! 
781   ENDIF 
782   !! 
783   !! skip OM_XXX diagnostics here 
784   !! 
785   IF( med_diag%TCO2%dgsave ) THEN 
786   CALL wrk_alloc( jpi, jpj, f_TDIC2d ) 
787   f_TDIC2d(:,:) = 0.0 !! 
788   ENDIF 
789   IF( med_diag%TALK%dgsave ) THEN 
790   CALL wrk_alloc( jpi, jpj, f_TALK2d ) 
791   f_TALK2d(:,:) = 0.0 !! 
792   ENDIF 
793   IF( med_diag%KW660%dgsave ) THEN 
794   CALL wrk_alloc( jpi, jpj, f_kw6602d ) 
795   f_kw6602d(:,:) = 0.0 !! 
796   ENDIF 
797   IF( med_diag%ATM_PP0%dgsave ) THEN 
798   CALL wrk_alloc( jpi, jpj, f_pp02d ) 
799   f_pp02d(:,:) = 0.0 !! 
800   ENDIF 
801   IF( med_diag%O2FLUX%dgsave ) THEN 
802   CALL wrk_alloc( jpi, jpj, f_o2flux2d ) 
803   f_o2flux2d(:,:) = 0.0 !! 
804   ENDIF 
805   IF( med_diag%O2SAT%dgsave ) THEN 
806   CALL wrk_alloc( jpi, jpj, f_o2sat2d ) 
807   f_o2sat2d(:,:) = 0.0 !! 
808   ENDIF 
809   !! 
810   !! skip XXX_CCD diagnostics here 
811   !! 
812   IF( med_diag%SFR_OCAL%dgsave ) THEN 
813   CALL wrk_alloc( jpi, jpj, sfr_ocal2d ) 
814   sfr_ocal2d(:,:) = 0.0 !! 
815   ENDIF 
816   IF( med_diag%SFR_OARG%dgsave ) THEN 
817   CALL wrk_alloc( jpi, jpj, sfr_oarg2d ) 
818   sfr_oarg2d(:,:) = 0.0 !! 
819   ENDIF 
820   !! 
821   !! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here 
822   !! 
823   IF( med_diag%IBEN_N%dgsave ) THEN 
824   CALL wrk_alloc( jpi, jpj, iben_n2d ) 
825   iben_n2d(:,:) = 0.0 !! 
826   ENDIF 
827   IF( med_diag%IBEN_FE%dgsave ) THEN 
828   CALL wrk_alloc( jpi, jpj, iben_fe2d ) 
829   iben_fe2d(:,:) = 0.0 !! 
830   ENDIF 
831   IF( med_diag%IBEN_C%dgsave ) THEN 
832   CALL wrk_alloc( jpi, jpj, iben_c2d ) 
833   iben_c2d(:,:) = 0.0 !! 
834   ENDIF 
835   IF( med_diag%IBEN_SI%dgsave ) THEN 
836   CALL wrk_alloc( jpi, jpj, iben_si2d ) 
837   iben_si2d(:,:) = 0.0 !! 
838   ENDIF 
839   IF( med_diag%IBEN_CA%dgsave ) THEN 
840   CALL wrk_alloc( jpi, jpj, iben_ca2d ) 
841   iben_ca2d(:,:) = 0.0 !! 
842   ENDIF 
843   IF( med_diag%OBEN_N%dgsave ) THEN 
844   CALL wrk_alloc( jpi, jpj, oben_n2d ) 
845   oben_n2d(:,:) = 0.0 !! 
846   ENDIF 
847   IF( med_diag%OBEN_FE%dgsave ) THEN 
848   CALL wrk_alloc( jpi, jpj, oben_fe2d ) 
849   oben_fe2d(:,:) = 0.0 !! 
850   ENDIF 
851   IF( med_diag%OBEN_C%dgsave ) THEN 
852   CALL wrk_alloc( jpi, jpj, oben_c2d ) 
853   oben_c2d(:,:) = 0.0 !! 
854   ENDIF 
855   IF( med_diag%OBEN_SI%dgsave ) THEN 
856   CALL wrk_alloc( jpi, jpj, oben_si2d ) 
857   oben_si2d(:,:) = 0.0 !! 
858   ENDIF 
859   IF( med_diag%OBEN_CA%dgsave ) THEN 
860   CALL wrk_alloc( jpi, jpj, oben_ca2d ) 
861   oben_ca2d(:,:) = 0.0 !! 
862   ENDIF 
863   !! 
864   !! skip BEN_XX diagnostics here 
865   !! 
866   IF( med_diag%RIV_N%dgsave ) THEN 
867   CALL wrk_alloc( jpi, jpj, rivn2d ) 
868   rivn2d(:,:) = 0.0 !! 
869   ENDIF 
870   IF( med_diag%RIV_SI%dgsave ) THEN 
871   CALL wrk_alloc( jpi, jpj, rivsi2d ) 
872   rivsi2d(:,:) = 0.0 !! 
873   ENDIF 
874   IF( med_diag%RIV_C%dgsave ) THEN 
875   CALL wrk_alloc( jpi, jpj, rivc2d ) 
876   rivc2d(:,:) = 0.0 !! 
877   ENDIF 
878   IF( med_diag%RIV_ALK%dgsave ) THEN 
879   CALL wrk_alloc( jpi, jpj, rivalk2d ) 
880   rivalk2d(:,:) = 0.0 !! 
881   ENDIF 
882   IF( med_diag%DETC%dgsave ) THEN 
883   CALL wrk_alloc( jpi, jpj, fslowc2d ) 
884   fslowc2d(:,:) = 0.0 !! 
885   ENDIF 
886   !! 
887   !! skip SDC_XXXX, INVTXXX diagnostics here 
888   !! 
889   IF( med_diag%LYSO_CA%dgsave ) THEN 
890   CALL wrk_alloc( jpi, jpj, lyso_ca2d ) 
891   lyso_ca2d(:,:) = 0.0 !! 
892   ENDIF 
893   !! 
894   !! skip COM_RESP diagnostic here 
895   !! 
896   IF( med_diag%PN_LLOSS%dgsave ) THEN 
897   CALL wrk_alloc( jpi, jpj, fdpn22d ) 
898   fdpn22d(:,:) = 0.0 !! 
899   ENDIF 
900   IF( med_diag%PD_LLOSS%dgsave ) THEN 
901   CALL wrk_alloc( jpi, jpj, fdpd22d ) 
902   fdpd22d(:,:) = 0.0 !! 
903   ENDIF 
904   IF( med_diag%ZI_LLOSS%dgsave ) THEN 
905   CALL wrk_alloc( jpi, jpj, fdzmi22d ) 
906   fdzmi22d(:,:) = 0.0 !! 
907   ENDIF 
908   IF( med_diag%ZE_LLOSS%dgsave ) THEN 
909   CALL wrk_alloc( jpi, jpj, fdzme22d ) 
910   fdzme22d(:,:) = 0.0 !! 
911   ENDIF 
912   IF( med_diag%ZI_MES_N%dgsave ) THEN 
913   CALL wrk_alloc( jpi, jpj, zimesn2d ) 
914   zimesn2d(:,:) = 0.0 !! 
915   ENDIF 
916   IF( med_diag%ZI_MES_D%dgsave ) THEN 
917   CALL wrk_alloc( jpi, jpj, zimesd2d ) 
918   zimesd2d(:,:) = 0.0 !! 
919   ENDIF 
920   IF( med_diag%ZI_MES_C%dgsave ) THEN 
921   CALL wrk_alloc( jpi, jpj, zimesc2d ) 
922   zimesc2d(:,:) = 0.0 !! 
923   ENDIF 
924   IF( med_diag%ZI_MESDC%dgsave ) THEN 
925   CALL wrk_alloc( jpi, jpj, zimesdc2d ) 
926   zimesdc2d(:,:) = 0.0 !! 
927   ENDIF 
928   IF( med_diag%ZI_EXCR%dgsave ) THEN 
929   CALL wrk_alloc( jpi, jpj, ziexcr2d ) 
930   ziexcr2d(:,:) = 0.0 !! 
931   ENDIF 
932   IF( med_diag%ZI_RESP%dgsave ) THEN 
933   CALL wrk_alloc( jpi, jpj, ziresp2d ) 
934   ziresp2d(:,:) = 0.0 !! 
935   ENDIF 
936   IF( med_diag%ZI_GROW%dgsave ) THEN 
937   CALL wrk_alloc( jpi, jpj, zigrow2d ) 
938   zigrow2d(:,:) = 0.0 !! 
939   ENDIF 
940   IF( med_diag%ZE_MES_N%dgsave ) THEN 
941   CALL wrk_alloc( jpi, jpj, zemesn2d ) 
942   zemesn2d(:,:) = 0.0 !! 
943   ENDIF 
944   IF( med_diag%ZE_MES_D%dgsave ) THEN 
945   CALL wrk_alloc( jpi, jpj, zemesd2d ) 
946   zemesd2d(:,:) = 0.0 !! 
947   ENDIF 
948   IF( med_diag%ZE_MES_C%dgsave ) THEN 
949   CALL wrk_alloc( jpi, jpj, zemesc2d ) 
950   zemesc2d(:,:) = 0.0 !! 
951   ENDIF 
952   IF( med_diag%ZE_MESDC%dgsave ) THEN 
953   CALL wrk_alloc( jpi, jpj, zemesdc2d ) 
954   zemesdc2d(:,:) = 0.0 !! 
955   ENDIF 
956   IF( med_diag%ZE_EXCR%dgsave ) THEN 
957   CALL wrk_alloc( jpi, jpj, zeexcr2d ) 
958   zeexcr2d(:,:) = 0.0 !! 
959   ENDIF 
960   IF( med_diag%ZE_RESP%dgsave ) THEN 
961   CALL wrk_alloc( jpi, jpj, zeresp2d ) 
962   zeresp2d(:,:) = 0.0 !! 
963   ENDIF 
964   IF( med_diag%ZE_GROW%dgsave ) THEN 
965   CALL wrk_alloc( jpi, jpj, zegrow2d ) 
966   zegrow2d(:,:) = 0.0 !! 
967   ENDIF 
968   IF( med_diag%MDETC%dgsave ) THEN 
969   CALL wrk_alloc( jpi, jpj, mdetc2d ) 
970   mdetc2d(:,:) = 0.0 !! 
971   ENDIF 
972   IF( med_diag%GMIDC%dgsave ) THEN 
973   CALL wrk_alloc( jpi, jpj, gmidc2d ) 
974   gmidc2d(:,:) = 0.0 !! 
975   ENDIF 
976   IF( med_diag%GMEDC%dgsave ) THEN 
977   CALL wrk_alloc( jpi, jpj, gmedc2d ) 
978   gmedc2d(:,:) = 0.0 !! 
979   ENDIF 
980   !! 
981   !! skip INT_XXX diagnostics here 
982   !! 
983   IF (jdms .eq. 1) THEN 
984   IF( med_diag%DMS_SURF%dgsave ) THEN 
985   CALL wrk_alloc( jpi, jpj, dms_surf2d ) 
986   dms_surf2d(:,:) = 0.0 !! 
987   ENDIF 
988   IF( med_diag%DMS_ANDR%dgsave ) THEN 
989   CALL wrk_alloc( jpi, jpj, dms_andr2d ) 
990   dms_andr2d(:,:) = 0.0 !! 
991   ENDIF 
992   IF( med_diag%DMS_SIMO%dgsave ) THEN 
993   CALL wrk_alloc( jpi, jpj, dms_simo2d ) 
994   dms_simo2d(:,:) = 0.0 !! 
995   ENDIF 
996   IF( med_diag%DMS_ARAN%dgsave ) THEN 
997   CALL wrk_alloc( jpi, jpj, dms_aran2d ) 
998   dms_aran2d(:,:) = 0.0 !! 
999   ENDIF 
1000   IF( med_diag%DMS_HALL%dgsave ) THEN 
1001   CALL wrk_alloc( jpi, jpj, dms_hall2d ) 
1002   dms_hall2d(:,:) = 0.0 !! 
1003   ENDIF 
1004   ENDIF 
1005   !! 
1006   !! AXY (24/11/16): extra MOCSY diagnostics, 2D 
1007   IF( med_diag%ATM_XCO2%dgsave ) THEN 
1008   CALL wrk_alloc( jpi, jpj, f_xco2a_2d ) 
1009   f_xco2a_2d(:,:) = 0.0 !! 
1010   ENDIF 
1011   IF( med_diag%OCN_FCO2%dgsave ) THEN 
1012   CALL wrk_alloc( jpi, jpj, f_fco2w_2d ) 
1013   f_fco2w_2d(:,:) = 0.0 !! 
1014   ENDIF 
1015   IF( med_diag%ATM_FCO2%dgsave ) THEN 
1016   CALL wrk_alloc( jpi, jpj, f_fco2a_2d ) 
1017   f_fco2a_2d(:,:) = 0.0 !! 
1018   ENDIF 
1019   IF( med_diag%OCN_RHOSW%dgsave ) THEN 
1020   CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d ) 
1021   f_ocnrhosw_2d(:,:) = 0.0 !! 
1022   ENDIF 
1023   IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
1024   CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d ) 
1025   f_ocnschco2_2d(:,:) = 0.0 !! 
1026   ENDIF 
1027   IF( med_diag%OCN_KWCO2%dgsave ) THEN 
1028   CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d ) 
1029   f_ocnkwco2_2d(:,:) = 0.0 !! 
1030   ENDIF 
1031   IF( med_diag%OCN_K0%dgsave ) THEN 
1032   CALL wrk_alloc( jpi, jpj, f_ocnk0_2d ) 
1033   f_ocnk0_2d(:,:) = 0.0 !! 
1034   ENDIF 
1035   IF( med_diag%CO2STARAIR%dgsave ) THEN 
1036   CALL wrk_alloc( jpi, jpj, f_co2starair_2d ) 
1037   f_co2starair_2d(:,:) = 0.0 !! 
1038   ENDIF 
1039   IF( med_diag%OCN_DPCO2%dgsave ) THEN 
1040   CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d ) 
1041   f_ocndpco2_2d(:,:) = 0.0 !! 
1042   ENDIF 
 738  !! 
 739  !! skip SEAFLRXX, MED_XXXX, INTFLX_XX, INT_XX, ML_XXX, OCAL_XXX, FE_XXXX, MED_XZE, WIND diagnostics here 
 740  !! 
 741  IF( med_diag%RR_0100%dgsave ) THEN 
 742  CALL wrk_alloc( jpi, jpj, ffastca2d ) 
 743  ffastca2d(:,:) = 0.0 !! 
 744  ENDIF 
 745  
 746  IF( med_diag%ATM_PCO2%dgsave ) THEN 
 747  CALL wrk_alloc( jpi, jpj, f_pco2a2d ) 
 748  f_pco2a2d(:,:) = 0.0 !! 
 749  ENDIF 
 750  !! 
 751  !! skip OCN_PH diagnostic here 
 752  !! 
 753  IF( med_diag%OCN_PCO2%dgsave ) THEN 
 754  CALL wrk_alloc( jpi, jpj, f_pco2w2d ) 
 755  f_pco2w2d(:,:) = 0.0 !! 
 756  ENDIF 
 757  !! 
 758  !! skip OCNH2CO3, OCN_HCO3, OCN_CO3 diagnostics here 
 759  !! 
 760  IF( med_diag%CO2FLUX%dgsave ) THEN 
 761  CALL wrk_alloc( jpi, jpj, f_co2flux2d ) 
 762  f_co2flux2d(:,:) = 0.0 !! 
 763  ENDIF 
 764  !! 
 765  !! skip OM_XXX diagnostics here 
 766  !! 
 767  IF( med_diag%TCO2%dgsave ) THEN 
 768  CALL wrk_alloc( jpi, jpj, f_TDIC2d ) 
 769  f_TDIC2d(:,:) = 0.0 !! 
 770  ENDIF 
 771  IF( med_diag%TALK%dgsave ) THEN 
 772  CALL wrk_alloc( jpi, jpj, f_TALK2d ) 
 773  f_TALK2d(:,:) = 0.0 !! 
 774  ENDIF 
 775  IF( med_diag%KW660%dgsave ) THEN 
 776  CALL wrk_alloc( jpi, jpj, f_kw6602d ) 
 777  f_kw6602d(:,:) = 0.0 !! 
 778  ENDIF 
 779  IF( med_diag%ATM_PP0%dgsave ) THEN 
 780  CALL wrk_alloc( jpi, jpj, f_pp02d ) 
 781  f_pp02d(:,:) = 0.0 !! 
 782  ENDIF 
 783  IF( med_diag%O2FLUX%dgsave ) THEN 
 784  CALL wrk_alloc( jpi, jpj, f_o2flux2d ) 
 785  f_o2flux2d(:,:) = 0.0 !! 
 786  ENDIF 
 787  IF( med_diag%O2SAT%dgsave ) THEN 
 788  CALL wrk_alloc( jpi, jpj, f_o2sat2d ) 
 789  f_o2sat2d(:,:) = 0.0 !! 
 790  ENDIF 
 791  !! 
 792  !! skip XXX_CCD diagnostics here 
 793  !! 
 794  IF( med_diag%SFR_OCAL%dgsave ) THEN 
 795  CALL wrk_alloc( jpi, jpj, sfr_ocal2d ) 
 796  sfr_ocal2d(:,:) = 0.0 !! 
 797  ENDIF 
 798  IF( med_diag%SFR_OARG%dgsave ) THEN 
 799  CALL wrk_alloc( jpi, jpj, sfr_oarg2d ) 
 800  sfr_oarg2d(:,:) = 0.0 !! 
 801  ENDIF 
 802  !! 
 803  !! skip XX_PROD, XX_CONS, O2_ANOX, RR_XXXX diagnostics here 
 804  !! 
 805  IF( med_diag%IBEN_N%dgsave ) THEN 
 806  CALL wrk_alloc( jpi, jpj, iben_n2d ) 
 807  iben_n2d(:,:) = 0.0 !! 
 808  ENDIF 
 809  IF( med_diag%IBEN_FE%dgsave ) THEN 
 810  CALL wrk_alloc( jpi, jpj, iben_fe2d ) 
 811  iben_fe2d(:,:) = 0.0 !! 
 812  ENDIF 
 813  IF( med_diag%IBEN_C%dgsave ) THEN 
 814  CALL wrk_alloc( jpi, jpj, iben_c2d ) 
 815  iben_c2d(:,:) = 0.0 !! 
 816  ENDIF 
 817  IF( med_diag%IBEN_SI%dgsave ) THEN 
 818  CALL wrk_alloc( jpi, jpj, iben_si2d ) 
 819  iben_si2d(:,:) = 0.0 !! 
 820  ENDIF 
 821  IF( med_diag%IBEN_CA%dgsave ) THEN 
 822  CALL wrk_alloc( jpi, jpj, iben_ca2d ) 
 823  iben_ca2d(:,:) = 0.0 !! 
 824  ENDIF 
 825  IF( med_diag%OBEN_N%dgsave ) THEN 
 826  CALL wrk_alloc( jpi, jpj, oben_n2d ) 
 827  oben_n2d(:,:) = 0.0 !! 
 828  ENDIF 
 829  IF( med_diag%OBEN_FE%dgsave ) THEN 
 830  CALL wrk_alloc( jpi, jpj, oben_fe2d ) 
 831  oben_fe2d(:,:) = 0.0 !! 
 832  ENDIF 
 833  IF( med_diag%OBEN_C%dgsave ) THEN 
 834  CALL wrk_alloc( jpi, jpj, oben_c2d ) 
 835  oben_c2d(:,:) = 0.0 !! 
 836  ENDIF 
 837  IF( med_diag%OBEN_SI%dgsave ) THEN 
 838  CALL wrk_alloc( jpi, jpj, oben_si2d ) 
 839  oben_si2d(:,:) = 0.0 !! 
 840  ENDIF 
 841  IF( med_diag%OBEN_CA%dgsave ) THEN 
 842  CALL wrk_alloc( jpi, jpj, oben_ca2d ) 
 843  oben_ca2d(:,:) = 0.0 !! 
 844  ENDIF 
 845  !! 
 846  !! skip BEN_XX diagnostics here 
 847  !! 
 848  IF( med_diag%RIV_N%dgsave ) THEN 
 849  CALL wrk_alloc( jpi, jpj, rivn2d ) 
 850  rivn2d(:,:) = 0.0 !! 
 851  ENDIF 
 852  IF( med_diag%RIV_SI%dgsave ) THEN 
 853  CALL wrk_alloc( jpi, jpj, rivsi2d ) 
 854  rivsi2d(:,:) = 0.0 !! 
 855  ENDIF 
 856  IF( med_diag%RIV_C%dgsave ) THEN 
 857  CALL wrk_alloc( jpi, jpj, rivc2d ) 
 858  rivc2d(:,:) = 0.0 !! 
 859  ENDIF 
 860  IF( med_diag%RIV_ALK%dgsave ) THEN 
 861  CALL wrk_alloc( jpi, jpj, rivalk2d ) 
 862  rivalk2d(:,:) = 0.0 !! 
 863  ENDIF 
 864  IF( med_diag%DETC%dgsave ) THEN 
 865  CALL wrk_alloc( jpi, jpj, fslowc2d ) 
 866  fslowc2d(:,:) = 0.0 !! 
 867  ENDIF 
 868  !! 
 869  !! skip SDC_XXXX, INVTXXX diagnostics here 
 870  !! 
 871  IF( med_diag%LYSO_CA%dgsave ) THEN 
 872  CALL wrk_alloc( jpi, jpj, lyso_ca2d ) 
 873  lyso_ca2d(:,:) = 0.0 !! 
 874  ENDIF 
 875  !! 
 876  !! skip COM_RESP diagnostic here 
 877  !! 
 878  IF( med_diag%PN_LLOSS%dgsave ) THEN 
 879  CALL wrk_alloc( jpi, jpj, fdpn22d ) 
 880  fdpn22d(:,:) = 0.0 !! 
 881  ENDIF 
 882  IF( med_diag%PD_LLOSS%dgsave ) THEN 
 883  CALL wrk_alloc( jpi, jpj, fdpd22d ) 
 884  fdpd22d(:,:) = 0.0 !! 
 885  ENDIF 
 886  IF( med_diag%ZI_LLOSS%dgsave ) THEN 
 887  CALL wrk_alloc( jpi, jpj, fdzmi22d ) 
 888  fdzmi22d(:,:) = 0.0 !! 
 889  ENDIF 
 890  IF( med_diag%ZE_LLOSS%dgsave ) THEN 
 891  CALL wrk_alloc( jpi, jpj, fdzme22d ) 
 892  fdzme22d(:,:) = 0.0 !! 
 893  ENDIF 
 894  IF( med_diag%ZI_MES_N%dgsave ) THEN 
 895  CALL wrk_alloc( jpi, jpj, zimesn2d ) 
 896  zimesn2d(:,:) = 0.0 !! 
 897  ENDIF 
 898  IF( med_diag%ZI_MES_D%dgsave ) THEN 
 899  CALL wrk_alloc( jpi, jpj, zimesd2d ) 
 900  zimesd2d(:,:) = 0.0 !! 
 901  ENDIF 
 902  IF( med_diag%ZI_MES_C%dgsave ) THEN 
 903  CALL wrk_alloc( jpi, jpj, zimesc2d ) 
 904  zimesc2d(:,:) = 0.0 !! 
 905  ENDIF 
 906  IF( med_diag%ZI_MESDC%dgsave ) THEN 
 907  CALL wrk_alloc( jpi, jpj, zimesdc2d ) 
 908  zimesdc2d(:,:) = 0.0 !! 
 909  ENDIF 
 910  IF( med_diag%ZI_EXCR%dgsave ) THEN 
 911  CALL wrk_alloc( jpi, jpj, ziexcr2d ) 
 912  ziexcr2d(:,:) = 0.0 !! 
 913  ENDIF 
 914  IF( med_diag%ZI_RESP%dgsave ) THEN 
 915  CALL wrk_alloc( jpi, jpj, ziresp2d ) 
 916  ziresp2d(:,:) = 0.0 !! 
 917  ENDIF 
 918  IF( med_diag%ZI_GROW%dgsave ) THEN 
 919  CALL wrk_alloc( jpi, jpj, zigrow2d ) 
 920  zigrow2d(:,:) = 0.0 !! 
 921  ENDIF 
 922  IF( med_diag%ZE_MES_N%dgsave ) THEN 
 923  CALL wrk_alloc( jpi, jpj, zemesn2d ) 
 924  zemesn2d(:,:) = 0.0 !! 
 925  ENDIF 
 926  IF( med_diag%ZE_MES_D%dgsave ) THEN 
 927  CALL wrk_alloc( jpi, jpj, zemesd2d ) 
 928  zemesd2d(:,:) = 0.0 !! 
 929  ENDIF 
 930  IF( med_diag%ZE_MES_C%dgsave ) THEN 
 931  CALL wrk_alloc( jpi, jpj, zemesc2d ) 
 932  zemesc2d(:,:) = 0.0 !! 
 933  ENDIF 
 934  IF( med_diag%ZE_MESDC%dgsave ) THEN 
 935  CALL wrk_alloc( jpi, jpj, zemesdc2d ) 
 936  zemesdc2d(:,:) = 0.0 !! 
 937  ENDIF 
 938  IF( med_diag%ZE_EXCR%dgsave ) THEN 
 939  CALL wrk_alloc( jpi, jpj, zeexcr2d ) 
 940  zeexcr2d(:,:) = 0.0 !! 
 941  ENDIF 
 942  IF( med_diag%ZE_RESP%dgsave ) THEN 
 943  CALL wrk_alloc( jpi, jpj, zeresp2d ) 
 944  zeresp2d(:,:) = 0.0 !! 
 945  ENDIF 
 946  IF( med_diag%ZE_GROW%dgsave ) THEN 
 947  CALL wrk_alloc( jpi, jpj, zegrow2d ) 
 948  zegrow2d(:,:) = 0.0 !! 
 949  ENDIF 
 950  IF( med_diag%MDETC%dgsave ) THEN 
 951  CALL wrk_alloc( jpi, jpj, mdetc2d ) 
 952  mdetc2d(:,:) = 0.0 !! 
 953  ENDIF 
 954  IF( med_diag%GMIDC%dgsave ) THEN 
 955  CALL wrk_alloc( jpi, jpj, gmidc2d ) 
 956  gmidc2d(:,:) = 0.0 !! 
 957  ENDIF 
 958  IF( med_diag%GMEDC%dgsave ) THEN 
 959  CALL wrk_alloc( jpi, jpj, gmedc2d ) 
 960  gmedc2d(:,:) = 0.0 !! 
 961  ENDIF 
 962  !! 
 963  !! skip INT_XXX diagnostics here 
 964  !! 
 965  IF (jdms .eq. 1) THEN 
 966  IF( med_diag%DMS_SURF%dgsave ) THEN 
 967  CALL wrk_alloc( jpi, jpj, dms_surf2d ) 
 968  dms_surf2d(:,:) = 0.0 !! 
 969  ENDIF 
 970  IF( med_diag%DMS_ANDR%dgsave ) THEN 
 971  CALL wrk_alloc( jpi, jpj, dms_andr2d ) 
 972  dms_andr2d(:,:) = 0.0 !! 
 973  ENDIF 
 974  IF( med_diag%DMS_SIMO%dgsave ) THEN 
 975  CALL wrk_alloc( jpi, jpj, dms_simo2d ) 
 976  dms_simo2d(:,:) = 0.0 !! 
 977  ENDIF 
 978  IF( med_diag%DMS_ARAN%dgsave ) THEN 
 979  CALL wrk_alloc( jpi, jpj, dms_aran2d ) 
 980  dms_aran2d(:,:) = 0.0 !! 
 981  ENDIF 
 982  IF( med_diag%DMS_HALL%dgsave ) THEN 
 983  CALL wrk_alloc( jpi, jpj, dms_hall2d ) 
 984  dms_hall2d(:,:) = 0.0 !! 
 985  ENDIF 
 986  ENDIF 
 987  !! 
 988  !! AXY (24/11/16): extra MOCSY diagnostics, 2D 
 989  IF( med_diag%ATM_XCO2%dgsave ) THEN 
 990  CALL wrk_alloc( jpi, jpj, f_xco2a_2d ) 
 991  f_xco2a_2d(:,:) = 0.0 !! 
 992  ENDIF 
 993  IF( med_diag%OCN_FCO2%dgsave ) THEN 
 994  CALL wrk_alloc( jpi, jpj, f_fco2w_2d ) 
 995  f_fco2w_2d(:,:) = 0.0 !! 
 996  ENDIF 
 997  IF( med_diag%ATM_FCO2%dgsave ) THEN 
 998  CALL wrk_alloc( jpi, jpj, f_fco2a_2d ) 
 999  f_fco2a_2d(:,:) = 0.0 !! 
 1000  ENDIF 
 1001  IF( med_diag%OCN_RHOSW%dgsave ) THEN 
 1002  CALL wrk_alloc( jpi, jpj, f_ocnrhosw_2d ) 
 1003  f_ocnrhosw_2d(:,:) = 0.0 !! 
 1004  ENDIF 
 1005  IF( med_diag%OCN_SCHCO2%dgsave ) THEN 
 1006  CALL wrk_alloc( jpi, jpj, f_ocnschco2_2d ) 
 1007  f_ocnschco2_2d(:,:) = 0.0 !! 
 1008  ENDIF 
 1009  IF( med_diag%OCN_KWCO2%dgsave ) THEN 
 1010  CALL wrk_alloc( jpi, jpj, f_ocnkwco2_2d ) 
 1011  f_ocnkwco2_2d(:,:) = 0.0 !! 
 1012  ENDIF 
 1013  IF( med_diag%OCN_K0%dgsave ) THEN 
 1014  CALL wrk_alloc( jpi, jpj, f_ocnk0_2d ) 
 1015  f_ocnk0_2d(:,:) = 0.0 !! 
 1016  ENDIF 
 1017  IF( med_diag%CO2STARAIR%dgsave ) THEN 
 1018  CALL wrk_alloc( jpi, jpj, f_co2starair_2d ) 
 1019  f_co2starair_2d(:,:) = 0.0 !! 
 1020  ENDIF 
 1021  IF( med_diag%OCN_DPCO2%dgsave ) THEN 
 1022  CALL wrk_alloc( jpi, jpj, f_ocndpco2_2d ) 
 1023  f_ocndpco2_2d(:,:) = 0.0 !! 
 1024  ENDIF 
1044   IF( med_diag%TPP3%dgsave ) THEN 
1045   CALL wrk_alloc( jpi, jpj, jpk, tpp3d ) 
1046   tpp3d(:,:,:) = 0.0 !! 
1047   ENDIF 
1048   IF( med_diag%DETFLUX3%dgsave ) THEN 
1049   CALL wrk_alloc( jpi, jpj, jpk, detflux3d ) 
1050   detflux3d(:,:,:) = 0.0 !! 
1051   ENDIF 
1052   IF( med_diag%REMIN3N%dgsave ) THEN 
1053   CALL wrk_alloc( jpi, jpj, jpk, remin3dn ) 
1054   remin3dn(:,:,:) = 0.0 !! 
1055   ENDIF 
1056   !! 
1057   !! AXY (10/11/16): CMIP6 diagnostics, 2D 
1058   !! JPALM  171116  put fgco2 alloc out of diag request 
1059   !! needed for coupling/passed through restart 
1060   !! IF( med_diag%FGCO2%dgsave ) THEN 
1061   CALL wrk_alloc( jpi, jpj, fgco2 ) 
1062   fgco2(:,:) = 0.0 !! 
1063   !! ENDIF 
1064   IF( med_diag%INTDISSIC%dgsave ) THEN 
1065   CALL wrk_alloc( jpi, jpj, intdissic ) 
1066   intdissic(:,:) = 0.0 !! 
1067   ENDIF 
1068   IF( med_diag%INTDISSIN%dgsave ) THEN 
1069   CALL wrk_alloc( jpi, jpj, intdissin ) 
1070   intdissin(:,:) = 0.0 !! 
1071   ENDIF 
1072   IF( med_diag%INTDISSISI%dgsave ) THEN 
1073   CALL wrk_alloc( jpi, jpj, intdissisi ) 
1074   intdissisi(:,:) = 0.0 !! 
1075   ENDIF 
1076   IF( med_diag%INTTALK%dgsave ) THEN 
1077   CALL wrk_alloc( jpi, jpj, inttalk ) 
1078   inttalk(:,:) = 0.0 !! 
1079   ENDIF 
1080   IF( med_diag%O2min%dgsave ) THEN 
1081   CALL wrk_alloc( jpi, jpj, o2min ) 
1082   o2min(:,:) = 1.e3 !! set to high value as we're looking for min(o2) 
1083   ENDIF 
1084   IF( med_diag%ZO2min%dgsave ) THEN 
1085   CALL wrk_alloc( jpi, jpj, zo2min ) 
1086   zo2min(:,:) = 0.0 !! 
1087   ENDIF 
1088   IF( med_diag%FBDDTALK%dgsave ) THEN 
1089   CALL wrk_alloc( jpi, jpj, fbddtalk ) 
1090   fbddtalk(:,:) = 0.0 !! 
1091   ENDIF 
1092   IF( med_diag%FBDDTDIC%dgsave ) THEN 
1093   CALL wrk_alloc( jpi, jpj, fbddtdic ) 
1094   fbddtdic(:,:) = 0.0 !! 
1095   ENDIF 
1096   IF( med_diag%FBDDTDIFE%dgsave ) THEN 
1097   CALL wrk_alloc( jpi, jpj, fbddtdife ) 
1098   fbddtdife(:,:) = 0.0 !! 
1099   ENDIF 
1100   IF( med_diag%FBDDTDIN%dgsave ) THEN 
1101   CALL wrk_alloc( jpi, jpj, fbddtdin ) 
1102   fbddtdin(:,:) = 0.0 !! 
1103   ENDIF 
1104   IF( med_diag%FBDDTDISI%dgsave ) THEN 
1105   CALL wrk_alloc( jpi, jpj, fbddtdisi ) 
1106   fbddtdisi(:,:) = 0.0 !! 
1107   ENDIF 
1108   !! 
1109   !! AXY (10/11/16): CMIP6 diagnostics, 3D 
1110   IF( med_diag%TPPD3%dgsave ) THEN 
1111   CALL wrk_alloc( jpi, jpj, jpk, tppd3 ) 
1112   tppd3(:,:,:) = 0.0 !! 
1113   ENDIF 
1114   IF( med_diag%BDDTALK3%dgsave ) THEN 
1115   CALL wrk_alloc( jpi, jpj, jpk, bddtalk3 ) 
1116   bddtalk3(:,:,:) = 0.0 !! 
1117   ENDIF 
1118   IF( med_diag%BDDTDIC3%dgsave ) THEN 
1119   CALL wrk_alloc( jpi, jpj, jpk, bddtdic3 ) 
1120   bddtdic3(:,:,:) = 0.0 !! 
1121   ENDIF 
1122   IF( med_diag%BDDTDIFE3%dgsave ) THEN 
1123   CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 ) 
1124   bddtdife3(:,:,:) = 0.0 !! 
1125   ENDIF 
1126   IF( med_diag%BDDTDIN3%dgsave ) THEN 
1127   CALL wrk_alloc( jpi, jpj, jpk, bddtdin3 ) 
1128   bddtdin3(:,:,:) = 0.0 !! 
1129   ENDIF 
1130   IF( med_diag%BDDTDISI3%dgsave ) THEN 
1131   CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 ) 
1132   bddtdisi3(:,:,:) = 0.0 !! 
1133   ENDIF 
1134   IF( med_diag%FD_NIT3%dgsave ) THEN 
1135   CALL wrk_alloc( jpi, jpj, jpk, fd_nit3 ) 
1136   fd_nit3(:,:,:) = 0.0 !! 
1137   ENDIF 
1138   IF( med_diag%FD_SIL3%dgsave ) THEN 
1139   CALL wrk_alloc( jpi, jpj, jpk, fd_sil3 ) 
1140   fd_sil3(:,:,:) = 0.0 !! 
1141   ENDIF 
1142   IF( med_diag%FD_CAR3%dgsave ) THEN 
1143   CALL wrk_alloc( jpi, jpj, jpk, fd_car3 ) 
1144   fd_car3(:,:,:) = 0.0 !! 
1145   ENDIF 
1146   IF( med_diag%FD_CAL3%dgsave ) THEN 
1147   CALL wrk_alloc( jpi, jpj, jpk, fd_cal3 ) 
1148   fd_cal3(:,:,:) = 0.0 !! 
1149   ENDIF 
1150   IF( med_diag%DCALC3%dgsave ) THEN 
1151   CALL wrk_alloc( jpi, jpj, jpk, dcalc3 ) 
1152   dcalc3(:,:,: ) = 0.0 !! 
1153   ENDIF 
1154   IF( med_diag%EXPC3%dgsave ) THEN 
1155   CALL wrk_alloc( jpi, jpj, jpk, expc3 ) 
1156   expc3(:,:,: ) = 0.0 !! 
1157   ENDIF 
1158   IF( med_diag%EXPN3%dgsave ) THEN 
1159   CALL wrk_alloc( jpi, jpj, jpk, expn3 ) 
1160   expn3(:,:,: ) = 0.0 !! 
1161   ENDIF 
1162   IF( med_diag%FEDISS3%dgsave ) THEN 
1163   CALL wrk_alloc( jpi, jpj, jpk, fediss3 ) 
1164   fediss3(:,:,: ) = 0.0 !! 
1165   ENDIF 
1166   IF( med_diag%FESCAV3%dgsave ) THEN 
1167   CALL wrk_alloc( jpi, jpj, jpk, fescav3 ) 
1168   fescav3(:,:,: ) = 0.0 !! 
1169   ENDIF 
1170   IF( med_diag%MIGRAZP3%dgsave ) THEN 
1171   CALL wrk_alloc( jpi, jpj, jpk, migrazp3 ) 
1172   migrazp3(:,:,: ) = 0.0 !! 
1173   ENDIF 
1174   IF( med_diag%MIGRAZD3%dgsave ) THEN 
1175   CALL wrk_alloc( jpi, jpj, jpk, migrazd3 ) 
1176   migrazd3(:,:,: ) = 0.0 !! 
1177   ENDIF 
1178   IF( med_diag%MEGRAZP3%dgsave ) THEN 
1179   CALL wrk_alloc( jpi, jpj, jpk, megrazp3 ) 
1180   megrazp3(:,:,: ) = 0.0 !! 
1181   ENDIF 
1182   IF( med_diag%MEGRAZD3%dgsave ) THEN 
1183   CALL wrk_alloc( jpi, jpj, jpk, megrazd3 ) 
1184   megrazd3(:,:,: ) = 0.0 !! 
1185   ENDIF 
1186   IF( med_diag%MEGRAZZ3%dgsave ) THEN 
1187   CALL wrk_alloc( jpi, jpj, jpk, megrazz3 ) 
1188   megrazz3(:,:,: ) = 0.0 !! 
1189   ENDIF 
1190   IF( med_diag%O2SAT3%dgsave ) THEN 
1191   CALL wrk_alloc( jpi, jpj, jpk, o2sat3 ) 
1192   o2sat3(:,:,: ) = 0.0 !! 
1193   ENDIF 
1194   IF( med_diag%PBSI3%dgsave ) THEN 
1195   CALL wrk_alloc( jpi, jpj, jpk, pbsi3 ) 
1196   pbsi3(:,:,: ) = 0.0 !! 
1197   ENDIF 
1198   IF( med_diag%PCAL3%dgsave ) THEN 
1199   CALL wrk_alloc( jpi, jpj, jpk, pcal3 ) 
1200   pcal3(:,:,: ) = 0.0 !! 
1201   ENDIF 
1202   IF( med_diag%REMOC3%dgsave ) THEN 
1203   CALL wrk_alloc( jpi, jpj, jpk, remoc3 ) 
1204   remoc3(:,:,: ) = 0.0 !! 
1205   ENDIF 
1206   IF( med_diag%PNLIMJ3%dgsave ) THEN 
1207   CALL wrk_alloc( jpi, jpj, jpk, pnlimj3 ) 
1208   pnlimj3(:,:,: ) = 0.0 !! 
1209   ENDIF 
1210   IF( med_diag%PNLIMN3%dgsave ) THEN 
1211   CALL wrk_alloc( jpi, jpj, jpk, pnlimn3 ) 
1212   pnlimn3(:,:,: ) = 0.0 !! 
1213   ENDIF 
1214   IF( med_diag%PNLIMFE3%dgsave ) THEN 
1215   CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3 ) 
1216   pnlimfe3(:,:,: ) = 0.0 !! 
1217   ENDIF 
1218   IF( med_diag%PDLIMJ3%dgsave ) THEN 
1219   CALL wrk_alloc( jpi, jpj, jpk, pdlimj3 ) 
1220   pdlimj3(:,:,: ) = 0.0 !! 
1221   ENDIF 
1222   IF( med_diag%PDLIMN3%dgsave ) THEN 
1223   CALL wrk_alloc( jpi, jpj, jpk, pdlimn3 ) 
1224   pdlimn3(:,:,: ) = 0.0 !! 
1225   ENDIF 
1226   IF( med_diag%PDLIMFE3%dgsave ) THEN 
1227   CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3 ) 
1228   pdlimfe3(:,:,: ) = 0.0 !! 
1229   ENDIF 
1230   IF( med_diag%PDLIMSI3%dgsave ) THEN 
1231   CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3 ) 
1232   pdlimsi3(:,:,: ) = 0.0 !! 
1233   ENDIF 
1234   
1235   ENDIF 
1236   !! lk_iomput 
1237   !! 
 1026  IF( med_diag%TPP3%dgsave ) THEN 
 1027  CALL wrk_alloc( jpi, jpj, jpk, tpp3d ) 
 1028  tpp3d(:,:,:) = 0.0 !! 
 1029  ENDIF 
 1030  IF( med_diag%DETFLUX3%dgsave ) THEN 
 1031  CALL wrk_alloc( jpi, jpj, jpk, detflux3d ) 
 1032  detflux3d(:,:,:) = 0.0 !! 
 1033  ENDIF 
 1034  IF( med_diag%REMIN3N%dgsave ) THEN 
 1035  CALL wrk_alloc( jpi, jpj, jpk, remin3dn ) 
 1036  remin3dn(:,:,:) = 0.0 !! 
 1037  ENDIF 
 1038  !! 
 1039  !! AXY (10/11/16): CMIP6 diagnostics, 2D 
 1040  !! JPALM  171116  put fgco2 alloc out of diag request 
 1041  !! needed for coupling/passed through restart 
 1042  !! IF( med_diag%FGCO2%dgsave ) THEN 
 1043  CALL wrk_alloc( jpi, jpj, fgco2 ) 
 1044  fgco2(:,:) = 0.0 !! 
 1045  !! ENDIF 
 1046  IF( med_diag%INTDISSIC%dgsave ) THEN 
 1047  CALL wrk_alloc( jpi, jpj, intdissic ) 
 1048  intdissic(:,:) = 0.0 !! 
 1049  ENDIF 
 1050  IF( med_diag%INTDISSIN%dgsave ) THEN 
 1051  CALL wrk_alloc( jpi, jpj, intdissin ) 
 1052  intdissin(:,:) = 0.0 !! 
 1053  ENDIF 
 1054  IF( med_diag%INTDISSISI%dgsave ) THEN 
 1055  CALL wrk_alloc( jpi, jpj, intdissisi ) 
 1056  intdissisi(:,:) = 0.0 !! 
 1057  ENDIF 
 1058  IF( med_diag%INTTALK%dgsave ) THEN 
 1059  CALL wrk_alloc( jpi, jpj, inttalk ) 
 1060  inttalk(:,:) = 0.0 !! 
 1061  ENDIF 
 1062  IF( med_diag%O2min%dgsave ) THEN 
 1063  CALL wrk_alloc( jpi, jpj, o2min ) 
 1064  o2min(:,:) = 1.e3 !! set to high value as we're looking for min(o2) 
 1065  ENDIF 
 1066  IF( med_diag%ZO2min%dgsave ) THEN 
 1067  CALL wrk_alloc( jpi, jpj, zo2min ) 
 1068  zo2min(:,:) = 0.0 !! 
 1069  ENDIF 
 1070  IF( med_diag%FBDDTALK%dgsave ) THEN 
 1071  CALL wrk_alloc( jpi, jpj, fbddtalk ) 
 1072  fbddtalk(:,:) = 0.0 !! 
 1073  ENDIF 
 1074  IF( med_diag%FBDDTDIC%dgsave ) THEN 
 1075  CALL wrk_alloc( jpi, jpj, fbddtdic ) 
 1076  fbddtdic(:,:) = 0.0 !! 
 1077  ENDIF 
 1078  IF( med_diag%FBDDTDIFE%dgsave ) THEN 
 1079  CALL wrk_alloc( jpi, jpj, fbddtdife ) 
 1080  fbddtdife(:,:) = 0.0 !! 
 1081  ENDIF 
 1082  IF( med_diag%FBDDTDIN%dgsave ) THEN 
 1083  CALL wrk_alloc( jpi, jpj, fbddtdin ) 
 1084  fbddtdin(:,:) = 0.0 !! 
 1085  ENDIF 
 1086  IF( med_diag%FBDDTDISI%dgsave ) THEN 
 1087  CALL wrk_alloc( jpi, jpj, fbddtdisi ) 
 1088  fbddtdisi(:,:) = 0.0 !! 
 1089  ENDIF 
 1090  !! 
 1091  !! AXY (10/11/16): CMIP6 diagnostics, 3D 
 1092  IF( med_diag%TPPD3%dgsave ) THEN 
 1093  CALL wrk_alloc( jpi, jpj, jpk, tppd3 ) 
 1094  tppd3(:,:,:) = 0.0 !! 
 1095  ENDIF 
 1096  IF( med_diag%BDDTALK3%dgsave ) THEN 
 1097  CALL wrk_alloc( jpi, jpj, jpk, bddtalk3 ) 
 1098  bddtalk3(:,:,:) = 0.0 !! 
 1099  ENDIF 
 1100  IF( med_diag%BDDTDIC3%dgsave ) THEN 
 1101  CALL wrk_alloc( jpi, jpj, jpk, bddtdic3 ) 
 1102  bddtdic3(:,:,:) = 0.0 !! 
 1103  ENDIF 
 1104  IF( med_diag%BDDTDIFE3%dgsave ) THEN 
 1105  CALL wrk_alloc( jpi, jpj, jpk, bddtdife3 ) 
 1106  bddtdife3(:,:,:) = 0.0 !! 
 1107  ENDIF 
 1108  IF( med_diag%BDDTDIN3%dgsave ) THEN 
 1109  CALL wrk_alloc( jpi, jpj, jpk, bddtdin3 ) 
 1110  bddtdin3(:,:,:) = 0.0 !! 
 1111  ENDIF 
 1112  IF( med_diag%BDDTDISI3%dgsave ) THEN 
 1113  CALL wrk_alloc( jpi, jpj, jpk, bddtdisi3 ) 
 1114  bddtdisi3(:,:,:) = 0.0 !! 
 1115  ENDIF 
 1116  IF( med_diag%FD_NIT3%dgsave ) THEN 
 1117  CALL wrk_alloc( jpi, jpj, jpk, fd_nit3 ) 
 1118  fd_nit3(:,:,:) = 0.0 !! 
 1119  ENDIF 
 1120  IF( med_diag%FD_SIL3%dgsave ) THEN 
 1121  CALL wrk_alloc( jpi, jpj, jpk, fd_sil3 ) 
 1122  fd_sil3(:,:,:) = 0.0 !! 
 1123  ENDIF 
 1124  IF( med_diag%FD_CAR3%dgsave ) THEN 
 1125  CALL wrk_alloc( jpi, jpj, jpk, fd_car3 ) 
 1126  fd_car3(:,:,:) = 0.0 !! 
 1127  ENDIF 
 1128  IF( med_diag%FD_CAL3%dgsave ) THEN 
 1129  CALL wrk_alloc( jpi, jpj, jpk, fd_cal3 ) 
 1130  fd_cal3(:,:,:) = 0.0 !! 
 1131  ENDIF 
 1132  IF( med_diag%DCALC3%dgsave ) THEN 
 1133  CALL wrk_alloc( jpi, jpj, jpk, dcalc3 ) 
 1134  dcalc3(:,:,: ) = 0.0 !! 
 1135  ENDIF 
 1136  IF( med_diag%EXPC3%dgsave ) THEN 
 1137  CALL wrk_alloc( jpi, jpj, jpk, expc3 ) 
 1138  expc3(:,:,: ) = 0.0 !! 
 1139  ENDIF 
 1140  IF( med_diag%EXPN3%dgsave ) THEN 
 1141  CALL wrk_alloc( jpi, jpj, jpk, expn3 ) 
 1142  expn3(:,:,: ) = 0.0 !! 
 1143  ENDIF 
 1144  IF( med_diag%FEDISS3%dgsave ) THEN 
 1145  CALL wrk_alloc( jpi, jpj, jpk, fediss3 ) 
 1146  fediss3(:,:,: ) = 0.0 !! 
 1147  ENDIF 
 1148  IF( med_diag%FESCAV3%dgsave ) THEN 
 1149  CALL wrk_alloc( jpi, jpj, jpk, fescav3 ) 
 1150  fescav3(:,:,: ) = 0.0 !! 
 1151  ENDIF 
 1152  IF( med_diag%MIGRAZP3%dgsave ) THEN 
 1153  CALL wrk_alloc( jpi, jpj, jpk, migrazp3 ) 
 1154  migrazp3(:,:,: ) = 0.0 !! 
 1155  ENDIF 
 1156  IF( med_diag%MIGRAZD3%dgsave ) THEN 
 1157  CALL wrk_alloc( jpi, jpj, jpk, migrazd3 ) 
 1158  migrazd3(:,:,: ) = 0.0 !! 
 1159  ENDIF 
 1160  IF( med_diag%MEGRAZP3%dgsave ) THEN 
 1161  CALL wrk_alloc( jpi, jpj, jpk, megrazp3 ) 
 1162  megrazp3(:,:,: ) = 0.0 !! 
 1163  ENDIF 
 1164  IF( med_diag%MEGRAZD3%dgsave ) THEN 
 1165  CALL wrk_alloc( jpi, jpj, jpk, megrazd3 ) 
 1166  megrazd3(:,:,: ) = 0.0 !! 
 1167  ENDIF 
 1168  IF( med_diag%MEGRAZZ3%dgsave ) THEN 
 1169  CALL wrk_alloc( jpi, jpj, jpk, megrazz3 ) 
 1170  megrazz3(:,:,: ) = 0.0 !! 
 1171  ENDIF 
 1172  IF( med_diag%O2SAT3%dgsave ) THEN 
 1173  CALL wrk_alloc( jpi, jpj, jpk, o2sat3 ) 
 1174  o2sat3(:,:,: ) = 0.0 !! 
 1175  ENDIF 
 1176  IF( med_diag%PBSI3%dgsave ) THEN 
 1177  CALL wrk_alloc( jpi, jpj, jpk, pbsi3 ) 
 1178  pbsi3(:,:,: ) = 0.0 !! 
 1179  ENDIF 
 1180  IF( med_diag%PCAL3%dgsave ) THEN 
 1181  CALL wrk_alloc( jpi, jpj, jpk, pcal3 ) 
 1182  pcal3(:,:,: ) = 0.0 !! 
 1183  ENDIF 
 1184  IF( med_diag%REMOC3%dgsave ) THEN 
 1185  CALL wrk_alloc( jpi, jpj, jpk, remoc3 ) 
 1186  remoc3(:,:,: ) = 0.0 !! 
 1187  ENDIF 
 1188  IF( med_diag%PNLIMJ3%dgsave ) THEN 
 1189  CALL wrk_alloc( jpi, jpj, jpk, pnlimj3 ) 
 1190  pnlimj3(:,:,: ) = 0.0 !! 
 1191  ENDIF 
 1192  IF( med_diag%PNLIMN3%dgsave ) THEN 
 1193  CALL wrk_alloc( jpi, jpj, jpk, pnlimn3 ) 
 1194  pnlimn3(:,:,: ) = 0.0 !! 
 1195  ENDIF 
 1196  IF( med_diag%PNLIMFE3%dgsave ) THEN 
 1197  CALL wrk_alloc( jpi, jpj, jpk, pnlimfe3 ) 
 1198  pnlimfe3(:,:,: ) = 0.0 !! 
 1199  ENDIF 
 1200  IF( med_diag%PDLIMJ3%dgsave ) THEN 
 1201  CALL wrk_alloc( jpi, jpj, jpk, pdlimj3 ) 
 1202  pdlimj3(:,:,: ) = 0.0 !! 
 1203  ENDIF 
 1204  IF( med_diag%PDLIMN3%dgsave ) THEN 
 1205  CALL wrk_alloc( jpi, jpj, jpk, pdlimn3 ) 
 1206  pdlimn3(:,:,: ) = 0.0 !! 
 1207  ENDIF 
 1208  IF( med_diag%PDLIMFE3%dgsave ) THEN 
 1209  CALL wrk_alloc( jpi, jpj, jpk, pdlimfe3 ) 
 1210  pdlimfe3(:,:,: ) = 0.0 !! 
 1211  ENDIF 
 1212  IF( med_diag%PDLIMSI3%dgsave ) THEN 
 1213  CALL wrk_alloc( jpi, jpj, jpk, pdlimsi3 ) 
 1214  pdlimsi3(:,:,: ) = 0.0 !! 
 1215  ENDIF 
 1216  ENDIF !! lk_iomput 
 1217  !! 
1263   CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' ) 
1264   endif 
 1233  enddo 
 1234  CALL ctl_stop( 'trcbio_medusa, NAN in incoming tracer field' ) 
 1235  endif 
 1236  ENDDO 
 1237  CALL flush(numout) 
 1238  # endif 
 1239  
 1240  # if defined key_debug_medusa 
 1241  IF ( lwp ) write (numout,*) 'trc_bio_medusa: variables initialised and checked' 
 1242  CALL flush(numout) 
 1243  # endif 
 1244  
 1245  # if defined key_roam 
 1246  !! 
 1247  !! calculate atmospheric pCO2 
 1248  !! 
 1249  !! 
 1250  # if defined key_axy_pi_co2 
 1251  f_xco2a = 284.725 !! OCMIP preindustrial pCO2 
 1252  # else 
 1253  f_xco2a = 284.725 !! OCMIP preindustrial pCO2 
 1254  # endif 
 1255  IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a 
 1256  # endif 
 1257  
 1258  # if defined key_debug_medusa 
 1259  IF ( lwp ) write (numout,*) 'trc_bio_medusa: ready for carbonate chemistry' 
 1260  IF ( lwp ) write (numout,*) 'trc_bio_medusa: kt = ', kt 
 1261  IF ( lwp ) write (numout,*) 'trc_bio_medusa: nittrc000 = ', nittrc000 
 1262  CALL flush(numout) 
 1263  # endif 
 1264  
 1265  !!====================================================================== 
 1266  !! AXY (07/04/17): possible subroutine block; ocean interior carbonate chemistry 
 1267  !!====================================================================== 
 1268  # if defined key_roam 
 1269  !! AXY (20/11/14): alter to call on first MEDUSA timestep and then every 
 1270  !! month (this is hardwired as 960 timesteps but should 
 1271  !! be calculated and done properly 
 1272  !! IF( kt == nit000 .or. mod(kt,1920) == 0 ) THEN 
 1273  !! IF( kt == nittrc000 .or. mod(kt,960) == 0 ) THEN 
 1274  !!============================= 
 1275  !! Jpalm  07102016  need to change carbchem frequency call : 
 1276  !! we don't want to call on the first timestep of all run submission, 
 1277  !! but only on the very first timestep, and then every month 
 1278  !! So we call on nittrc000 if not restarted run, 
 1279  !! else if one month after last call. 
 1280  !! assume one month is 30d > 3600*24*30 : 2592000s 
 1281  !! try to call carbchem at 1st month's tmstp : x * 30d + 1*rdt(i.e: mod = rdt) 
 1282  !! ++ need to pass carbchem output var through restarts 
 1283  IF ( ( kt == nittrc000 .AND. .NOT.ln_rsttr ) .OR. mod(kt*rdt,2592000.) == rdt ) THEN 
 1284  !! 
 1285  !! Calculate the carbonate chemistry for the whole ocean on the first 
 1286  !! simulation timestep and every month subsequently; the resulting 3D 
 1287  !! field of omega calcite is used to determine the depth of the CCD 
 1288  !! 
 1289  !! 
 1290  IF(lwp) WRITE(numout,*) ' MEDUSA calculating all carbonate chemistry at kt =', kt 
 1291  CALL flush(numout) 
 1292  !! blank flags 
 1293  i2_omcal(:,:) = 0 
 1294  i2_omarg(:,:) = 0 
 1295  !! loop over 3D space 
 1296  DO jk = 1,jpk 
 1297  DO jj = 2,jpjm1 
 1298  DO ji = 2,jpim1 
 1299  !! OPEN wet point IF..THEN loop 
 1300  if (tmask(ji,jj,jk).eq.1) then 
 1301  IF ( lk_oasis ) THEN 
 1302  f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 
 1303  ENDIF 
 1304  !! AXY (06/04/17): where am I? 
 1305  flatx = gphit(ji,jj) 
 1306  !! do carbonate chemistry 
 1307  !! 
 1308  fdep2 = fsdept(ji,jj,jk) !! set up level midpoint 
 1309  !! AXY (28/11/16): seafloor depth; previously mbathy(ji,jj)  1, now mbathy(ji,jj) 
 1310  jmbathy = mbathy(ji,jj) 
 1311  !! 
 1312  !! set up required state variables 
 1313  zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
 1314  zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 
 1315  ztmp = tsn(ji,jj,jk,jp_tem) !! temperature 
 1316  zsal = tsn(ji,jj,jk,jp_sal) !! salinity 
 1317  zsil = max(0.,trn(ji,jj,jk,jpsil)) !! silicic acid 
 1318  zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
 1319  !! 
 1320  !! AXY (28/02/14): check input fields 
 1321  if (ztmp .lt. 3.0 .or. ztmp .gt. 40.0 ) then 
 1322  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 3D, ', & 
 1323  tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 
 1324  ji, ',', jj, ',', jk, ') at time', kt 
 1325  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 3D, ', & 
 1326  tsn(ji,jj,jk,jp_tem), ' > ', tsb(ji,jj,jk,jp_tem) 
 1327  ztmp = tsb(ji,jj,jk,jp_tem) !! temperature 
 1328  endif 
 1329  if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then 
 1330  IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 3D, ', & 
 1331  tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & 
 1332  ji, ',', jj, ',', jk, ') at time', kt 
 1333  endif 
 1334  !! 
 1335  !! blank input variables not used at this stage (they relate to airsea flux) 
 1336  f_kw660 = 1.0 
 1337  f_pp0 = 1.0 
 1338  !! 
 1339  !! calculate carbonate chemistry at grid cell midpoint 
 1340  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY2 carbonate 
 1341  !! chemistry package 
 1342  CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 
 1343  f_pp0, fdep2, flatx, f_kw660, f_xco2a, 1, & ! inputs 
 1344  f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 
 1345  f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 
 1346  f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 
 1347  f_co2starair, f_co2flux, f_dpco2 ) ! outputs 
 1348  !! 
 1349  f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 > umol / kg 
 1350  f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 > ueq / kg 
 1351  f_dcf = f_rhosw 
 1352  !! 
 1353  !! store 3D outputs 
 1354  f3_pH(ji,jj,jk) = f_ph 
 1355  f3_h2co3(ji,jj,jk) = f_h2co3 
 1356  f3_hco3(ji,jj,jk) = f_hco3 
 1357  f3_co3(ji,jj,jk) = f_co3 
 1358  f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 
 1359  f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 
 1360  !! 
 1361  !! CCD calculation: calcite 
 1362  if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then 
 1363  if (jk .eq. 1) then 
 1364  f2_ccd_cal(ji,jj) = fdep2 
 1365  else 
 1366  fq0 = f3_omcal(ji,jj,jk1)  f_omcal(ji,jj) 
 1367  fq1 = f3_omcal(ji,jj,jk1)  1.0 
 1368  fq2 = fq1 / (fq0 + tiny(fq0)) 
 1369  fq3 = fdep2  fsdept(ji,jj,jk1) 
 1370  fq4 = fq2 * fq3 
 1371  f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk1) + fq4 
 1372  endif 
 1373  i2_omcal(ji,jj) = 1 
 1374  endif 
 1375  if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
 1376  !! reached seafloor and still no dissolution; set to seafloor (Wpoint) 
 1377  f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1) 
 1378  i2_omcal(ji,jj) = 1 
 1379  endif 
 1380  !! 
 1381  !! CCD calculation: aragonite 
 1382  if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then 
 1383  if (jk .eq. 1) then 
 1384  f2_ccd_arg(ji,jj) = fdep2 
 1385  else 
 1386  fq0 = f3_omarg(ji,jj,jk1)  f_omarg(ji,jj) 
 1387  fq1 = f3_omarg(ji,jj,jk1)  1.0 
 1388  fq2 = fq1 / (fq0 + tiny(fq0)) 
 1389  fq3 = fdep2  fsdept(ji,jj,jk1) 
 1390  fq4 = fq2 * fq3 
 1391  f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk1) + fq4 
 1392  endif 
 1393  i2_omarg(ji,jj) = 1 
 1394  endif 
 1395  if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
 1396  !! reached seafloor and still no dissolution; set to seafloor (Wpoint) 
 1397  f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1) 
 1398  i2_omarg(ji,jj) = 1 
 1399  endif 
 1400  endif 
 1401  ENDDO 
 1402  ENDDO 
1389   endif 
1390   !! 
1391   !! blank input variables not used at this stage (they relate to airsea flux) 
1392   f_kw660 = 1.0 
1393   f_pp0 = 1.0 
1394   !! 
1395   !! calculate carbonate chemistry at grid cell midpoint 
1396   # if defined key_mocsy 
1397   !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY2 carbonate 
1398   !! chemistry package 
1399   CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 
1400   f_pp0, fdep2, gphit(ji,jj), f_kw660, f_xco2a, 1, & ! inputs 
1401   f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 
1402   f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 
1403   f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 
1404   f_co2starair, f_co2flux, f_dpco2 ) ! outputs 
1405   !! 
1406   f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 > umol / kg 
1407   f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 > ueq / kg 
1408   f_dcf = f_rhosw 
1409   # else 
1410   !! AXY (22/06/15): use old PML carbonate chemistry package (the 
1411   !! MEDUSA2 default) 
1412   CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, fdep2, f_kw660, & ! inputs 
1413   f_xco2a, f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 
1414   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters) ! outputs 
1415   !! 
1416   !! AXY (28/02/14): check output fields 
1417   if (iters .eq. 25) then 
1418   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: 3D ITERS WARNING, ', & 
1419   iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
1420   endif 
1421   # endif 
1422   !! 
1423   !! store 3D outputs 
1424   f3_pH(ji,jj,jk) = f_ph 
1425   f3_h2co3(ji,jj,jk) = f_h2co3 
1426   f3_hco3(ji,jj,jk) = f_hco3 
1427   f3_co3(ji,jj,jk) = f_co3 
1428   f3_omcal(ji,jj,jk) = f_omcal(ji,jj) 
1429   f3_omarg(ji,jj,jk) = f_omarg(ji,jj) 
1430   !! 
1431   !! CCD calculation: calcite 
1432   if (i2_omcal(ji,jj) .eq. 0 .and. f_omcal(ji,jj) .lt. 1.0) then 
1433   if (jk .eq. 1) then 
1434   f2_ccd_cal(ji,jj) = fdep2 
1435   else 
1436   fq0 = f3_omcal(ji,jj,jk1)  f_omcal(ji,jj) 
1437   fq1 = f3_omcal(ji,jj,jk1)  1.0 
1438   fq2 = fq1 / (fq0 + tiny(fq0)) 
1439   fq3 = fdep2  fsdept(ji,jj,jk1) 
1440   fq4 = fq2 * fq3 
1441   f2_ccd_cal(ji,jj) = fsdept(ji,jj,jk1) + fq4 
1442   endif 
1443   i2_omcal(ji,jj) = 1 
1444   endif 
1445   if ( i2_omcal(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
1446   !! reached seafloor and still no dissolution; set to seafloor (Wpoint) 
1447   f2_ccd_cal(ji,jj) = fsdepw(ji,jj,jk+1) 
1448   i2_omcal(ji,jj) = 1 
1449   endif 
1450   !! 
1451   !! CCD calculation: aragonite 
1452   if (i2_omarg(ji,jj) .eq. 0 .and. f_omarg(ji,jj) .lt. 1.0) then 
1453   if (jk .eq. 1) then 
1454   f2_ccd_arg(ji,jj) = fdep2 
1455   else 
1456   fq0 = f3_omarg(ji,jj,jk1)  f_omarg(ji,jj) 
1457   fq1 = f3_omarg(ji,jj,jk1)  1.0 
1458   fq2 = fq1 / (fq0 + tiny(fq0)) 
1459   fq3 = fdep2  fsdept(ji,jj,jk1) 
1460   fq4 = fq2 * fq3 
1461   f2_ccd_arg(ji,jj) = fsdept(ji,jj,jk1) + fq4 
1462   endif 
1463   i2_omarg(ji,jj) = 1 
1464   endif 
1465   if ( i2_omarg(ji,jj) .eq. 0 .and. jk .eq. jmbathy ) then 
1466   !! reached seafloor and still no dissolution; set to seafloor (Wpoint) 
1467   f2_ccd_arg(ji,jj) = fsdepw(ji,jj,jk+1) 
1468   i2_omarg(ji,jj) = 1 
1469   endif 
1470   endif 
1471   ENDDO 
1472   ENDDO 
1473   ENDDO 
1474   ENDIF 
1475   # endif 
1476   
1477   # if defined key_debug_medusa 
1478   IF (lwp) write (numout,*) 'trc_bio_medusa: ready for full domain calculations' 
1479   CALL flush(numout) 
1480   # endif 
1481   
1482   !! 
1483   !! MEDUSA has unified equation through the water column 
1484   !! (Diff. from LOBSTER which has two sets: bio and nonbio layers) 
1485   !! Statement below in LOBSTER is different: DO jk = 1, jpkbm1 
1486   !! 
1487   !! 
1488   !! NOTE: the ordering of the loops below differs from that of some other 
1489   !! models; looping over the vertical dimension is the outermost loop and 
1490   !! this complicates some calculations (e.g. storage of vertical fluxes 
1491   !! that can otherwise be done via a singular variable require 2D fields 
1492   !! here); however, these issues are relatively easily resolved, but the 
1493   !! loops CANNOT be reordered without potentially causing code efficiency 
1494   !! problems (e.g. array indexing means that reordering the loops would 
1495   !! require skipping between widelyspaced memory location; potentially 
1496   !! outside those immediately cached) 
1497   !! 
1498   !! OPEN vertical loop 
1499   DO jk = 1,jpk 
1500   !! OPEN horizontal loops 
1501   DO jj = 2,jpjm1 
1502   DO ji = 2,jpim1 
1503   !! OPEN wet point IF..THEN loop 
1504   if (tmask(ji,jj,jk).eq.1) then 
1505   !!====================================================================== 
1506   !! SETUP LOCAL GRID CELL 
1507   !!====================================================================== 
1508   !! 
1509   !! 
1510   !! Some notes on grid vertical structure 
1511   !!  fsdepw(ji,jj,jk) is the depth of the upper surface of level jk 
1512   !!  fsde3w(ji,jj,jk) is *approximately* the midpoint of level jk 
1513   !!  fse3t(ji,jj,jk) is the thickness of level jk 
1514   !! 
1515   !! 
1516   !! AXY (11/12/08): set up level thickness 
1517   fthk = fse3t(ji,jj,jk) 
1518   !! AXY (25/02/10): set up level depth (top of level) 
1519   fdep = fsdepw(ji,jj,jk) 
1520   !! AXY (01/03/10): set up level depth (bottom of level) 
1521   fdep1 = fdep + fthk 
1522   !! AXY (28/11/16): local seafloor depth 
1523   !! previously mbathy(ji,jj)  1, now mbathy(ji,jj) 
1524   jmbathy = mbathy(ji,jj) 
1525   !! 
1526   !! set up model tracers 
1527   !! negative values of state variables are not allowed to 
1528   !! contribute to the calculated fluxes 
1529   zchn = max(0.,trn(ji,jj,jk,jpchn)) !! nondiatom chlorophyll 
1530   zchd = max(0.,trn(ji,jj,jk,jpchd)) !! diatom chlorophyll 
1531   zphn = max(0.,trn(ji,jj,jk,jpphn)) !! nondiatoms 
1532   zphd = max(0.,trn(ji,jj,jk,jpphd)) !! diatoms 
1533   zpds = max(0.,trn(ji,jj,jk,jppds)) !! diatom silicon 
1534   !! AXY (28/01/10): probably need to take account of chl/biomass connection 
1535   if (zchn.eq.0.) zphn = 0. 
1536   if (zchd.eq.0.) zphd = 0. 
1537   if (zphn.eq.0.) zchn = 0. 
1538   if (zphd.eq.0.) zchd = 0. 
1539   !! AXY (23/01/14): duh  why did I forget diatom silicon? 
1540   if (zpds.eq.0.) zphd = 0. 
1541   if (zphd.eq.0.) zpds = 0. 
1542   zzmi = max(0.,trn(ji,jj,jk,jpzmi)) !! microzooplankton 
1543   zzme = max(0.,trn(ji,jj,jk,jpzme)) !! mesozooplankton 
1544   zdet = max(0.,trn(ji,jj,jk,jpdet)) !! detrital nitrogen 
1545   zdin = max(0.,trn(ji,jj,jk,jpdin)) !! dissolved inorganic nitrogen 
1546   zsil = max(0.,trn(ji,jj,jk,jpsil)) !! dissolved silicic acid 
1547   zfer = max(0.,trn(ji,jj,jk,jpfer)) !! dissolved "iron" 
1548   # if defined key_roam 
1549   zdtc = max(0.,trn(ji,jj,jk,jpdtc)) !! detrital carbon 
1550   zdic = max(0.,trn(ji,jj,jk,jpdic)) !! dissolved inorganic carbon 
1551   zalk = max(0.,trn(ji,jj,jk,jpalk)) !! alkalinity 
1552   zoxy = max(0.,trn(ji,jj,jk,jpoxy)) !! oxygen 
1553   # if defined key_axy_carbchem && defined key_mocsy 
1554   zpho = max(0.,trn(ji,jj,jk,jpdin)) / 16.0 !! phosphate via DIN and Redfield 
1555   # endif 
1556   !! 
1557   !! also need physical parameters for gas exchange calculations 
1558   ztmp = tsn(ji,jj,jk,jp_tem) 
1559   zsal = tsn(ji,jj,jk,jp_sal) 
1560   !! 
1561   !! AXY (28/02/14): check input fields 
1562   if (ztmp .lt. 3.0 .or. ztmp .gt. 40.0 ) then 
1563   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T WARNING 2D, ', & 
1564   tsb(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_tem), ' at (', & 
1565   ji, ',', jj, ',', jk, ') at time', kt 
1566   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: T SWITCHING 2D, ', & 
1567   tsn(ji,jj,jk,jp_tem), ' > ', tsb(ji,jj,jk,jp_tem) 
1568   ztmp = tsb(ji,jj,jk,jp_tem) !! temperature 
1569   endif 
1570   if (zsal .lt. 0.0 .or. zsal .gt. 45.0 ) then 
1571   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: S WARNING 2D, ', & 
1572   tsb(ji,jj,jk,jp_sal), tsn(ji,jj,jk,jp_sal), ' at (', & 
1573   ji, ',', jj, ',', jk, ') at time', kt 
1574   endif 
 1507  endif 
1610   # if defined key_debug_medusa 
1611   !! report state variable values 
1612   if (idf.eq.1.AND.idfval.eq.1) then 
1613   IF (lwp) write (numout,*) '' 
1614   IF (lwp) write (numout,*) 'fthk(',jk,') = ', fthk 
1615   IF (lwp) write (numout,*) 'zphn(',jk,') = ', zphn 
1616   IF (lwp) write (numout,*) 'zphd(',jk,') = ', zphd 
1617   IF (lwp) write (numout,*) 'zpds(',jk,') = ', zpds 
1618   IF (lwp) write (numout,*) 'zzmi(',jk,') = ', zzmi 
1619   IF (lwp) write (numout,*) 'zzme(',jk,') = ', zzme 
1620   IF (lwp) write (numout,*) 'zdet(',jk,') = ', zdet 
1621   IF (lwp) write (numout,*) 'zdin(',jk,') = ', zdin 
1622   IF (lwp) write (numout,*) 'zsil(',jk,') = ', zsil 
1623   IF (lwp) write (numout,*) 'zfer(',jk,') = ', zfer 
1624   # if defined key_roam 
1625   IF (lwp) write (numout,*) 'zdtc(',jk,') = ', zdtc 
1626   IF (lwp) write (numout,*) 'zdic(',jk,') = ', zdic 
1627   IF (lwp) write (numout,*) 'zalk(',jk,') = ', zalk 
1628   IF (lwp) write (numout,*) 'zoxy(',jk,') = ', zoxy 
 1543  !! sum tracers for inventory checks 
 1544  IF( lk_iomput ) THEN 
 1545  IF ( med_diag%INVTN%dgsave ) THEN 
 1546  ftot_n(ji,jj) = ftot_n(ji,jj) + & 
 1547  (fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) ) 
 1548  ENDIF 
 1549  IF ( med_diag%INVTSI%dgsave ) THEN 
 1550  ftot_si(ji,jj) = ftot_si(ji,jj) + & 
 1551  (fthk * ( zpds + zsil ) ) 
 1552  ENDIF 
 1553  IF ( med_diag%INVTFE%dgsave ) THEN 
 1554  ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 
 1555  (fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) ) 
 1556  ENDIF 
 1557  # if defined key_roam 
 1558  IF ( med_diag%INVTC%dgsave ) THEN 
 1559  ftot_c(ji,jj) = ftot_c(ji,jj) + & 
 1560  (fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + & 
 1561  (xthetazmi * zzmi) + (xthetazme * zzme) + zdtc + & 
 1562  zdic ) ) 
 1563  ENDIF 
 1564  IF ( med_diag%INVTALK%dgsave ) THEN 
 1565  ftot_a(ji,jj) = ftot_a(ji,jj) + (fthk * ( zalk ) ) 
 1566  ENDIF 
 1567  IF ( med_diag%INVTO2%dgsave ) THEN 
 1568  ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fthk * ( zoxy ) ) 
 1569  ENDIF 
 1570  !! 
 1571  !! AXY (10/11/16): CMIP6 diagnostics 
 1572  IF ( med_diag%INTDISSIC%dgsave ) THEN 
 1573  intdissic(ji,jj) = intdissic(ji,jj) + (fthk * zdic) 
 1574  ENDIF 
 1575  IF ( med_diag%INTDISSIN%dgsave ) THEN 
 1576  intdissin(ji,jj) = intdissin(ji,jj) + (fthk * zdin) 
 1577  ENDIF 
 1578  IF ( med_diag%INTDISSISI%dgsave ) THEN 
 1579  intdissisi(ji,jj) = intdissisi(ji,jj) + (fthk * zsil) 
 1580  ENDIF 
 1581  IF ( med_diag%INTTALK%dgsave ) THEN 
 1582  inttalk(ji,jj) = inttalk(ji,jj) + (fthk * zalk) 
 1583  ENDIF 
 1584  IF ( med_diag%O2min%dgsave ) THEN 
 1585  if ( zoxy < o2min(ji,jj) ) then 
 1586  o2min(ji,jj) = zoxy 
 1587  IF ( med_diag%ZO2min%dgsave ) THEN 
 1588  zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint 
 1589  ENDIF 
 1590  endif 
 1591  ENDIF 
 1592  # endif 
 1593  ENDIF 
 1594  
 1595  CALL flush(numout) 
 1596  
 1597  !!====================================================================== 
 1598  !! LOCAL GRID CELL CALCULATIONS 
 1599  !!====================================================================== 
 1600  !! 
 1601  !!====================================================================== 
 1602  !! AXY (07/04/17): possible subroutine block; airsea gas exchange 
 1603  !!====================================================================== 
 1604  # if defined key_roam 
 1605  if ( jk .eq. 1 ) then 
 1606  !! 
 1607  !! Airsea gas exchange 
 1608  !! 
 1609  !! 
 1610  !! AXY (17/07/14): zwind_i and zwind_j do not exist in this 
 1611  !! version of NEMO because it does not include 
 1612  !! the SBC changes that our local version has 
 1613  !! for accessing the HadGEM2 forcing; they 
 1614  !! could be added, but an alternative approach 
 1615  !! is to make use of wndm from oce_trc.F90 
 1616  !! which is wind speed at 10m (which is what 
 1617  !! is required here; this may need to be 
 1618  !! revisited when MEDUSA properly interacts 
 1619  !! with UKESM1 physics 
 1620  !! 
 1621  f_wind = wndm(ji,jj) 
 1622  IF ( lk_oasis ) THEN 
 1623  f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 
 1624  ENDIF 
 1625  !! 
 1626  !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
 1627  !! in MEDUSA, the gas transfer velocity used in the carbon 
 1628  !! and oxygen cycles has been harmonised and is calculated 
 1629  !! by the same function here; this harmonisation includes 
 1630  !! changes to the PML carbonate chemistry scheme so that 
 1631  !! it too makes use of the same gas transfer velocity; the 
 1632  !! preferred parameterisation of this is Wanninkhof (2014), 
 1633  !! option 7 
 1634  !! 
 1635  # if defined key_debug_medusa 
 1636  IF ( lwp ) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
 1637  CALL flush(numout) 
 1638  # endif 
 1639  CALL gas_transfer( f_wind, 1, 7, & ! inputs 
 1640  f_kw660 ) ! outputs 
 1641  # if defined key_debug_medusa 
 1642  IF ( lwp ) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
 1643  CALL flush(numout) 
 1644  # endif 
 1645  !! 
 1646  !! air pressure (atm); ultimately this will use air pressure at the base 
 1647  !! of the UKESM1 atmosphere 
 1648  !! 
 1649  f_pp0 = 1.0 
 1650  !! 
 1651  !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp 
 1652  !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
 1653  !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 
 1654  !! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind =', f_wind 
 1655  !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) 
 1656  !! 
 1657  # if defined key_axy_carbchem 
 1658  !! 
 1659  !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY2 carbonate 
 1660  !! chemistry package; note that depth is set to 
 1661  !! zero in this call 
 1662  CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 
 1663  f_pp0, 0.0, flatx, f_kw660, f_xco2a, 1, & ! inputs 
 1664  f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 
 1665  f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 
 1666  f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 
 1667  f_co2starair, f_co2flux, f_dpco2 ) ! outputs 
 1668  !! 
 1669  f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 > umol / kg 
 1670  f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 > ueq / kg 
 1671  f_dcf = f_rhosw 
 1672  # else 
 1673  !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
 1674  !! quasisensible alternatives 
 1675  f_ph = 8.1 
 1676  f_pco2w = f_xco2a 
 1677  f_h2co3 = 0.005 * zdic 
 1678  f_hco3 = 0.865 * zdic 
 1679  f_co3 = 0.130 * zdic 
 1680  f_omcal(ji,jj) = 4. 
 1681  f_omarg(ji,jj) = 2. 
 1682  f_co2flux = 0. 
 1683  f_TDIC = zdic 
 1684  f_TALK = zalk 
 1685  f_dcf = 1.026 
 1686  f_henry = 1. 
 1687  !! AXY (23/06/15): add in some extra MOCSY diagnostics 
 1688  f_fco2w = f_xco2a 
 1689  f_BetaD = 1. 
 1690  f_rhosw = 1.026 
 1691  f_opres = 0. 
 1692  f_insitut = ztmp 
 1693  f_pco2atm = f_xco2a 
 1694  f_fco2atm = f_xco2a 
 1695  f_schmidtco2 = 660. 
 1696  f_kwco2 = 0. 
 1697  f_K0 = 0. 
 1698  f_co2starair = f_xco2a 
 1699  f_dpco2 = 0. 
1630   endif 
1631   # endif 
1632   
1633   # if defined key_debug_medusa 
1634   if (idf.eq.1.AND.idfval.eq.1.AND.jk.eq.1) then 
1635   IF (lwp) write (numout,*) '' 
1636   IF (lwp) write (numout,*) 'dust = ', dust(ji,jj) 
1637   endif 
1638   # endif 
1639   
1640   !! sum tracers for inventory checks 
1641   IF( lk_iomput ) THEN 
1642   IF ( med_diag%INVTN%dgsave ) THEN 
1643   ftot_n(ji,jj) = ftot_n(ji,jj) + & 
1644   (fthk * ( zphn + zphd + zzmi + zzme + zdet + zdin ) ) 
1645   ENDIF 
1646   IF ( med_diag%INVTSI%dgsave ) THEN 
1647   ftot_si(ji,jj) = ftot_si(ji,jj) + & 
1648   (fthk * ( zpds + zsil ) ) 
1649   ENDIF 
1650   IF ( med_diag%INVTFE%dgsave ) THEN 
1651   ftot_fe(ji,jj) = ftot_fe(ji,jj) + & 
1652   (fthk * ( xrfn * ( zphn + zphd + zzmi + zzme + zdet ) + zfer ) ) 
1653   ENDIF 
1654   # if defined key_roam 
1655   IF ( med_diag%INVTC%dgsave ) THEN 
1656   ftot_c(ji,jj) = ftot_c(ji,jj) + & 
1657   (fthk * ( (xthetapn * zphn) + (xthetapd * zphd) + & 
1658   (xthetazmi * zzmi) + (xthetazme * zzme) + zdtc + & 
1659   zdic ) ) 
1660   ENDIF 
1661   IF ( med_diag%INVTALK%dgsave ) THEN 
1662   ftot_a(ji,jj) = ftot_a(ji,jj) + (fthk * ( zalk ) ) 
1663   ENDIF 
1664   IF ( med_diag%INVTO2%dgsave ) THEN 
1665   ftot_o2(ji,jj) = ftot_o2(ji,jj) + (fthk * ( zoxy ) ) 
1666   ENDIF 
1667   !! 
1668   !! AXY (10/11/16): CMIP6 diagnostics 
1669   IF ( med_diag%INTDISSIC%dgsave ) THEN 
1670   intdissic(ji,jj) = intdissic(ji,jj) + (fthk * zdic) 
1671   ENDIF 
1672   IF ( med_diag%INTDISSIN%dgsave ) THEN 
1673   intdissin(ji,jj) = intdissin(ji,jj) + (fthk * zdin) 
1674   ENDIF 
1675   IF ( med_diag%INTDISSISI%dgsave ) THEN 
1676   intdissisi(ji,jj) = intdissisi(ji,jj) + (fthk * zsil) 
1677   ENDIF 
1678   IF ( med_diag%INTTALK%dgsave ) THEN 
1679   inttalk(ji,jj) = inttalk(ji,jj) + (fthk * zalk) 
1680   ENDIF 
1681   IF ( med_diag%O2min%dgsave ) THEN 
1682   if ( zoxy < o2min(ji,jj) ) then 
1683   o2min(ji,jj) = zoxy 
1684   IF ( med_diag%ZO2min%dgsave ) THEN 
1685   zo2min(ji,jj) = (fdep + fdep1) / 2. !! layer midpoint 
1686   ENDIF 
1687   endif 
1688   ENDIF 
1689   # endif 
1690   ENDIF 
1691   
1692   CALL flush(numout) 
1693   
1694   !!====================================================================== 
1695   !! LOCAL GRID CELL CALCULATIONS 
1696   !!====================================================================== 
1697   !! 
1698   # if defined key_roam 
1699   if ( jk .eq. 1 ) then 
1700   !! 
1701   !! Airsea gas exchange 
1702   !! 
1703   !! 
1704   !! AXY (17/07/14): zwind_i and zwind_j do not exist in this 
1705   !! version of NEMO because it does not include 
1706   !! the SBC changes that our local version has 
1707   !! for accessing the HadGEM2 forcing; they 
1708   !! could be added, but an alternative approach 
1709   !! is to make use of wndm from oce_trc.F90 
1710   !! which is wind speed at 10m (which is what 
1711   !! is required here; this may need to be 
1712   !! revisited when MEDUSA properly interacts 
1713   !! with UKESM1 physics 
1714   !! 
1715   f_wind = wndm(ji,jj) 
1716   IF (lk_oasis) THEN 
1717   f_xco2a = PCO2a_in_cpl(ji,jj) !! use 2D atm xCO2 from atm coupling 
1718   ENDIF 
1719   !! 
1720   !! AXY (23/06/15): as part of an effort to update the carbonate chemistry 
1721   !! in MEDUSA, the gas transfer velocity used in the carbon 
1722   !! and oxygen cycles has been harmonised and is calculated 
1723   !! by the same function here; this harmonisation includes 
1724   !! changes to the PML carbonate chemistry scheme so that 
1725   !! it too makes use of the same gas transfer velocity; the 
1726   !! preferred parameterisation of this is Wanninkhof (2014), 
1727   !! option 7 
1728   !! 
 1701  !! 
 1702  !! mmol/m2/s > mmol/m3/d; correct for seaice; divide through by layer thickness 
 1703  f_co2flux = (1.  fr_i(ji,jj)) * f_co2flux * 86400. / fthk 
 1704  !! 
 1705  !! oxygen (O2); OCMIP2 code 
 1706  !! AXY (23/06/15): amend input list for oxygen to account for common gas 
 1707  !! transfer velocity 
 1708  !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs 
 1709  !! f_kw660, f_o2flux, f_o2sat ) ! outputs 
 1710  CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs 
 1711  f_kwo2, f_o2flux, f_o2sat ) ! outputs 
 1712  !! 
 1713  !! mmol/m2/s > mol/m3/d; correct for seaice; divide through by layer thickness 
 1714  f_o2flux = (1.  fr_i(ji,jj)) * f_o2flux * 86400. / fthk 
 1715  !! 
 1716  !! Jpalm (082014) 
 1717  !! DMS surface concentration calculation; initialy added using METOFFICE subroutine 
 1718  !! airsea flux calculated in atmospheric chemistry from atm and ocn concentrations 
 1719  !! 
 1720  !! AXY (13/03/15): this is amended to calculate all of the DMS 
 1721  !! estimates examined during UKESM1 (see comments 
 1722  !! in trcdms_medusa.F90) 
 1723  !! 
 1724  !! AXY (28/03/17): amended to pass DIN limitation instead of DIN concentration; 
 1725  !! accounts for differences in nutrient halfsaturations; changes 
 1726  !! also made in trc_dms_medusa 
 1727  !! 
 1728  IF (jdms .eq. 1) THEN 
 1729  !! 
 1730  !! calculate weighted halfsaturation for DIN uptake 
 1731  dms_wtkn = ((zphn * xnln) + (zphd * xnld)) / (zphn + zphd) 
 1732  !! 
 1733  !! feed in correct inputs 
 1734  if (jdms_input .eq. 0) then 
 1735  !! use instantaneous inputs 
 1736  dms_nlim = zdin / (zdin + dms_wtkn) 
 1737  !! 
 1738  CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), dms_nlim, & ! inputs 
 1739  dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 
 1740  else 
 1741  !! use dielaverage inputs 
 1742  dms_nlim = zn_dms_din(ji,jj) / (zn_dms_din(ji,jj) + dms_wtkn) 
 1743  !! 
 1744  CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs 
 1745  zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), dms_nlim, & ! inputs 
 1746  dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 
 1747  endif 
 1748  !! 
 1749  !! assign correct output to variable passed to atmosphere 
 1750  if (jdms_model .eq. 1) then 
 1751  dms_surf = dms_andr 
 1752  elseif (jdms_model .eq. 2) then 
 1753  dms_surf = dms_simo 
 1754  elseif (jdms_model .eq. 3) then 
 1755  dms_surf = dms_aran 
 1756  elseif (jdms_model .eq. 4) then 
 1757  dms_surf = dms_hall 
 1758  endif 
 1759  !! 
 1760  !! 2D diag through iom_use 
 1761  IF( lk_iomput ) THEN 
 1762  IF( med_diag%DMS_SURF%dgsave ) THEN 
 1763  dms_surf2d(ji,jj) = dms_surf 
 1764  ENDIF 
 1765  IF( med_diag%DMS_ANDR%dgsave ) THEN 
 1766  dms_andr2d(ji,jj) = dms_andr 
 1767  ENDIF 
 1768  IF( med_diag%DMS_SIMO%dgsave ) THEN 
 1769  dms_simo2d(ji,jj) = dms_simo 
 1770  ENDIF 
 1771  IF( med_diag%DMS_ARAN%dgsave ) THEN 
 1772  dms_aran2d(ji,jj) = dms_aran 
 1773  ENDIF 
 1774  IF( med_diag%DMS_HALL%dgsave ) THEN 
 1775  dms_hall2d(ji,jj) = dms_hall 
 1776  ENDIF 
1730   IF (lwp) write (numout,*) 'trc_bio_medusa: entering gas_transfer' 
1731   CALL flush(numout) 
1732   # endif 
1733   CALL gas_transfer( f_wind, 1, 7, & ! inputs 
1734   f_kw660 ) ! outputs 
1735   # if defined key_debug_medusa 
1736   IF (lwp) write (numout,*) 'trc_bio_medusa: exiting gas_transfer' 
1737   CALL flush(numout) 
1738   # endif 
1739   !! 
1740   !! air pressure (atm); ultimately this will use air pressure at the base 
1741   !! of the UKESM1 atmosphere 
1742   !! 
1743   f_pp0 = 1.0 
1744   !! 
1745   !! IF(lwp) WRITE(numout,*) ' MEDUSA ztmp =', ztmp 
1746   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_i =', zwind_i(ji,jj) 
1747   !! IF(lwp) WRITE(numout,*) ' MEDUSA zwind_j =', zwind_j(ji,jj) 
1748   !! IF(lwp) WRITE(numout,*) ' MEDUSA f_wind =', f_wind 
1749   !! IF(lwp) WRITE(numout,*) ' MEDUSA fr_i =', fr_i(ji,jj) 
1750   !! 
1751   # if defined key_axy_carbchem 
1752   # if defined key_mocsy 
1753   !! 
1754   !! AXY (22/06/15): use Orr & Epitalon (2015) MOCSY2 carbonate 
1755   !! chemistry package; note that depth is set to 
1756   !! zero in this call 
1757   CALL mocsy_interface( ztmp, zsal, zalk, zdic, zsil, zpho, & ! inputs 
1758   f_pp0, 0.0, gphit(ji,jj), f_kw660, f_xco2a, 1, & ! inputs 
1759   f_ph, f_pco2w, f_fco2w, f_h2co3, f_hco3, f_co3, f_omarg(ji,jj), & ! outputs 
1760   f_omcal(ji,jj), f_BetaD, f_rhosw, f_opres, f_insitut, & ! outputs 
1761   f_pco2atm, f_fco2atm, f_schmidtco2, f_kwco2, f_K0, & ! outputs 
1762   f_co2starair, f_co2flux, f_dpco2 ) ! outputs 
1763   !! 
1764   f_TDIC = (zdic / f_rhosw) * 1000. ! mmol / m3 > umol / kg 
1765   f_TALK = (zalk / f_rhosw) * 1000. ! meq / m3 > ueq / kg 
1766   f_dcf = f_rhosw 
1767   # else 
1768   iters = 0 
1769   !! 
1770   !! carbon dioxide (CO2); Jerry Blackford code (ostensibly OCMIP2, but not) 
1771   CALL trc_co2_medusa( ztmp, zsal, zdic, zalk, 0.0, f_kw660, f_xco2a, & ! inputs 
1772   f_ph, f_pco2w, f_h2co3, f_hco3, f_co3, f_omcal(ji,jj), & ! outputs 
1773   f_omarg(ji,jj), f_co2flux, f_TDIC, f_TALK, f_dcf, f_henry, iters ) ! outputs 
1774   !! 
1775   !! AXY (09/01/14): removed iteration and NaN checks; these have 
1776   !! been moved to trc_co2_medusa together with a 
1777   !! fudge that amends erroneous values (this is 
1778   !! intended to be a temporary fudge!); the 
1779   !! output warnings are retained here so that 
1780   !! failure position can be determined 
1781   if (iters .eq. 25) then 
1782   IF(lwp) WRITE(numout,*) ' trc_bio_medusa: ITERS WARNING, ', & 
1783   iters, ' AT (', ji, ', ', jj, ', ', jk, ') AT ', kt 
1784   endif 
1785   # endif 
1786   # else 
1787   !! AXY (18/04/13): switch off carbonate chemistry calculations; provide 
1788   !! quasisensible alternatives 
1789   f_ph = 8.1 
1790   f_pco2w = f_xco2a 
1791   f_h2co3 = 0.005 * zdic 
1792   f_hco3 = 0.865 * zdic 
1793   f_co3 = 0.130 * zdic 
1794   f_omcal(ji,jj) = 4. 
1795   f_omarg(ji,jj) = 2. 
1796   f_co2flux = 0. 
1797   f_TDIC = zdic 
1798   f_TALK = zalk 
1799   f_dcf = 1.026 
1800   f_henry = 1. 
1801   !! AXY (23/06/15): add in some extra MOCSY diagnostics 
1802   f_fco2w = f_xco2a 
1803   f_BetaD = 1. 
1804   f_rhosw = 1.026 
1805   f_opres = 0. 
1806   f_insitut = ztmp 
1807   f_pco2atm = f_xco2a 
1808   f_fco2atm = f_xco2a 
1809   f_schmidtco2 = 660. 
1810   f_kwco2 = 0. 
1811   f_K0 = 0. 
1812   f_co2starair = f_xco2a 
1813   f_dpco2 = 0. 
1814   # endif 
1815   !! 
1816   !! mmol/m2/s > mmol/m3/d; correct for seaice; divide through by layer thickness 
1817   f_co2flux = (1.  fr_i(ji,jj)) * f_co2flux * 86400. / fthk 
1818   !! 
1819   !! oxygen (O2); OCMIP2 code 
1820   !! AXY (23/06/15): amend input list for oxygen to account for common gas 
1821   !! transfer velocity 
1822   !! CALL trc_oxy_medusa( ztmp, zsal, f_uwind, f_vwind, f_pp0, zoxy / 1000., fthk, & ! inputs 
1823   !! f_kw660, f_o2flux, f_o2sat ) ! outputs 
1824   CALL trc_oxy_medusa( ztmp, zsal, f_kw660, f_pp0, zoxy, & ! inputs 
1825   f_kwo2, f_o2flux, f_o2sat ) ! outputs 
1826   !! 
1827   !! mmol/m2/s > mol/m3/d; correct for seaice; divide through by layer thickness 
1828   f_o2flux = (1.  fr_i(ji,jj)) * f_o2flux * 86400. / fthk 
1829   !! 
1830   !! Jpalm (082014) 
1831   !! DMS surface concentration calculation 
1832   !! initialy added for UKESM1 model. 
1833   !! using METOFFICE subroutine. 
1834   !! DMS module only needs Chl concentration and MLD 
1835   !! to get an aproximate value of DMS concentration. 
1836   !! airsea fluxes are calculated by atmospheric chemitry model 
1837   !! from atm and ocsurface concentrations. 
1838   !! 
1839   !! AXY (13/03/15): this is amended to calculate all of the DMS 
1840   !! estimates examined during UKESM1 (see comments 
1841   !! in trcdms_medusa.F90) 
1842   !! 
1843   IF (jdms .eq. 1) THEN 
1844   !! 
1845   !! feed in correct inputs 
1846   if (jdms_input .eq. 0) then 
1847   !! use instantaneous inputs 
1848   CALL trc_dms_medusa( zchn, zchd, hmld(ji,jj), qsr(ji,jj), zdin, & ! inputs 
1849   dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 
1850   else 
1851   !! use dielaverage inputs 
1852   CALL trc_dms_medusa( zn_dms_chn(ji,jj), zn_dms_chd(ji,jj), & ! inputs 
1853   zn_dms_mld(ji,jj), zn_dms_qsr(ji,jj), zn_dms_din(ji,jj), & ! inputs 
1854   dms_andr, dms_simo, dms_aran, dms_hall ) ! outputs 
1855   endif 
1856   !! 
1857   !! assign correct output to variable passed to atmosphere 
1858   if (jdms_model .eq. 1) then 
1859   dms_surf = dms_andr 
1860   elseif (jdms_model .eq. 2) then 
1861   dms_surf = dms_simo 
1862   elseif (jdms_model .eq. 3) then 
1863   dms_surf = dms_aran 
1864   elseif (jdms_model .eq. 4) then 
1865   dms_surf = dms_hall 
1866   endif 
1867   !! 
1868   !! 2D diag through iom_use 
1869   IF( lk_iomput ) THEN 
1870   IF( med_diag%DMS_SURF%dgsave ) THEN 
1871   dms_surf2d(ji,jj) = dms_surf 
1872   ENDIF 
1873   IF( med_diag%DMS_ANDR%dgsave ) THEN 
1874   dms_andr2d(ji,jj) = dms_andr 
1875   ENDIF 
1876   IF( med_diag%DMS_SIMO%dgsave ) THEN 
1877   dms_simo2d(ji,jj) = dms_simo 
1878   ENDIF 
1879   IF( med_diag%DMS_ARAN%dgsave ) THEN 
1880   dms_aran2d(ji,jj) = dms_aran 
1881   ENDIF 
1882   IF( med_diag%DMS_HALL%dgsave ) THEN 
1883   dms_hall2d(ji,jj) = dms_hall 
1884   ENDIF 
1885   # if defined key_debug_medusa 
1886   IF (lwp) write (numout,*) 'trc_bio_medusa: finish calculating dms' 
1887   CALL flush(numout) 
 1778  IF ( lwp ) write (numout,*) 'trc_bio_medusa: finish calculating dms' 
 1779  CALL flush(numout) 
1971   if ( jk .eq. 1 ) then 
1972   !! 
1973   !! River inputs 
1974   !! 
1975   !! 
1976   !! runoff comes in as kg / m2 / s 
1977   !! used and written out as m3 / m2 / d (= m / d) 
1978   !! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 
1979   !! 
1980   !! AXY (17/07/14): the compiler doesn't like this line for some reason; 
1981   !! as MEDUSA doesn't even use runoff for riverine inputs, 
1982   !! a temporary solution is to switch off runoff entirely 
1983   !! here; again, this change is one of several that will 
1984   !! need revisiting once MEDUSA has bedded down in UKESM1; 
1985   !! particularly so if the land scheme provides information 
1986   !! concerning nutrient fluxes 
1987   !! 
1988   !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 
1989   f_runoff(ji,jj) = 0.0 
1990   !! 
1991   !! nutrients are added via rivers to the model in one of two ways: 
1992   !! 1. via river concentration; i.e. the average nutrient concentration 
1993   !! of a river water is described by a spatial file, and this is 
1994   !! multiplied by runoff to give a nutrient flux 
1995   !! 2. via direct river flux; i.e. the average nutrient flux due to 
1996   !! rivers is described by a spatial file, and this is simply applied 
1997   !! as a direct nutrient flux (i.e. it does not relate or respond to 
1998   !! model runoff) 
1999   !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 
2000   !! alkalinity are derived from continentscale DIC estimates (Huang et al., 
2001   !! 2012) and some Arctic river alkalinity estimates (Katya?) 
2002   !! 
2003   !! as of 19/07/12, riverine nutrients can now be spread vertically across 
2004   !! several grid cells rather than just poured into the surface box; this 
2005   !! block of code is still executed, however, to set up the total amounts 
2006   !! of nutrient entering via rivers 
2007   !! 
2008   !! nitrogen 
2009   if (jriver_n .eq. 1) then 
2010   !! river concentration specified; use runoff to calculate input 
2011   f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 
2012   elseif (jriver_n .eq. 2) then 
2013   !! river flux specified; independent of runoff 
2014   f_riv_n(ji,jj) = riv_n(ji,jj) 
2015   endif 
2016   !! 
2017   !! silicon 
2018   if (jriver_si .eq. 1) then 
2019   !! river concentration specified; use runoff to calculate input 
2020   f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 
2021   elseif (jriver_si .eq. 2) then 
2022   !! river flux specified; independent of runoff 
2023   f_riv_si(ji,jj) = riv_si(ji,jj) 
2024   endif 
2025   !! 
2026   !! carbon 
2027   if (jriver_c .eq. 1) then 
2028   !! river concentration specified; use runoff to calculate input 
2029   f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 
2030   elseif (jriver_c .eq. 2) then 
2031   !! river flux specified; independent of runoff 
2032   f_riv_c(ji,jj) = riv_c(ji,jj) 
2033   endif 
2034   !! 
2035   !! alkalinity 
2036   if (jriver_alk .eq. 1) then 
2037   !! river concentration specified; use runoff to calculate input 
2038   f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 
2039   elseif (jriver_alk .eq. 2) then 
2040   !! river flux specified; independent of runoff 
2041   f_riv_alk(ji,jj) = riv_alk(ji,jj) 
2042   endif 
2043   
2044   endif 
2045   
2046   !! 
2047   !! Chlorophyll calculations 
2048   !! 
2049   !! 
2050   !! nondiatoms 
2051   if (zphn.GT.rsmall) then 
2052   fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn))) 
2053   faln = xaln * fthetan 
2054   else 
2055   fthetan = 0. 
2056   faln = 0. 
2057   endif 
2058   !! 
2059   !! diatoms 
2060   if (zphd.GT.rsmall) then 
2061   fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd))) 
2062   fald = xald * fthetad 
2063   else 
2064   fthetad = 0. 
2065   fald = 0. 
2066   endif 
 1863  !!====================================================================== 
 1864  !! AXY (07/04/17): possible subroutine block; riverine inputs (or delete; it's unused presently) 
 1865  !!====================================================================== 
 1866  if ( jk .eq. 1 ) then 
 1867  !! 
 1868  !! River inputs 
 1869  !! 
 1870  !! 
 1871  !! runoff comes in as kg / m2 / s 
 1872  !! used and written out as m3 / m2 / d (= m / d) 
 1873  !! where 1000 kg / m2 / d = 1 m3 / m2 / d = 1 m / d 
 1874  !! 
 1875  !! AXY (17/07/14): the compiler doesn't like this line for some reason; 
 1876  !! as MEDUSA doesn't even use runoff for riverine inputs, 
 1877  !! a temporary solution is to switch off runoff entirely 
 1878  !! here; again, this change is one of several that will 
 1879  !! need revisiting once MEDUSA has bedded down in UKESM1; 
 1880  !! particularly so if the land scheme provides information 
 1881  !! concerning nutrient fluxes 
 1882  !! 
 1883  !! f_runoff(ji,jj) = sf_rnf(1)%fnow(ji,jj,1) / 1000. * 60. * 60. * 24. 
 1884  f_runoff(ji,jj) = 0.0 
 1885  !! 
 1886  !! nutrients are added via rivers to the model in one of two ways: 
 1887  !! 1. via river concentration; i.e. the average nutrient concentration 
 1888  !! of a river water is described by a spatial file, and this is 
 1889  !! multiplied by runoff to give a nutrient flux 
 1890  !! 2. via direct river flux; i.e. the average nutrient flux due to 
 1891  !! rivers is described by a spatial file, and this is simply applied 
 1892  !! as a direct nutrient flux (i.e. it does not relate or respond to 
 1893  !! model runoff) 
 1894  !! nutrient fields are derived from the GlobalNEWS 2 database; carbon and 
 1895  !! alkalinity are derived from continentscale DIC estimates (Huang et al., 
 1896  !! 2012) and some Arctic river alkalinity estimates (Katya?) 
 1897  !! 
 1898  !! as of 19/07/12, riverine nutrients can now be spread vertically across 
 1899  !! several grid cells rather than just poured into the surface box; this 
 1900  !! block of code is still executed, however, to set up the total amounts 
 1901  !! of nutrient entering via rivers 
 1902  !! 
 1903  !! nitrogen 
 1904  if (jriver_n .eq. 1) then 
 1905  !! river concentration specified; use runoff to calculate input 
 1906  f_riv_n(ji,jj) = f_runoff(ji,jj) * riv_n(ji,jj) 
 1907  elseif (jriver_n .eq. 2) then 
 1908  !! river flux specified; independent of runoff 
 1909  f_riv_n(ji,jj) = riv_n(ji,jj) 
 1910  endif 
 1911  !! 
 1912  !! silicon 
 1913  if (jriver_si .eq. 1) then 
 1914  !! river concentration specified; use runoff to calculate input 
 1915  f_riv_si(ji,jj) = f_runoff(ji,jj) * riv_si(ji,jj) 
 1916  elseif (jriver_si .eq. 2) then 
 1917  !! river flux specified; independent of runoff 
 1918  f_riv_si(ji,jj) = riv_si(ji,jj) 
 1919  endif 
 1920  !! 
 1921  !! carbon 
 1922  if (jriver_c .eq. 1) then 
 1923  !! river concentration specified; use runoff to calculate input 
 1924  f_riv_c(ji,jj) = f_runoff(ji,jj) * riv_c(ji,jj) 
 1925  elseif (jriver_c .eq. 2) then 
 1926  !! river flux specified; independent of runoff 
 1927  f_riv_c(ji,jj) = riv_c(ji,jj) 
 1928  endif 
 1929  !! 
 1930  !! alkalinity 
 1931  if (jriver_alk .eq. 1) then 
 1932  !! river concentration specified; use runoff to calculate input 
 1933  f_riv_alk(ji,jj) = f_runoff(ji,jj) * riv_alk(ji,jj) 
 1934  elseif (jriver_alk .eq. 2) then 
 1935  !! river flux specified; independent of runoff 
 1936  f_riv_alk(ji,jj) = riv_alk(ji,jj) 
 1937  endif 
 1938  
 1939  endif 
 1940  
 1941  !!====================================================================== 
 1942  !! AXY (07/04/17): possible subroutine block; phytoplankton growth 
 1943  !!====================================================================== 
 1944  
 1945  !! 
 1946  !! Chlorophyll calculations 
 1947  !! 
 1948  !! 
 1949  !! nondiatoms 
 1950  if (zphn.GT.rsmall) then 
 1951  fthetan = max(tiny(zchn), (zchn * xxi) / (zphn + tiny(zphn))) 
 1952  faln = xaln * fthetan 
 1953  else 
 1954  fthetan = 0. 
 1955  faln = 0. 
 1956  endif 
 1957  !! 
 1958  !! diatoms 
 1959  if (zphd.GT.rsmall) then 
 1960  fthetad = max(tiny(zchd), (zchd * xxi) / (zphd + tiny(zphd))) 
 1961  fald = xald * fthetad 
 1962  else 
 1963  fthetad = 0. 
 1964  fald = 0. 
 1965  endif 
 1966  
 1967  !! 
 1968  !! Phytoplankton light limitation 
 1969  !! 
 1970  !! 
 1971  !! It is assumed xpar is the depthaveraged (vertical layer) PAR 
 1972  !! Light limitation (check selfshading) in W/m2 
 1973  !! 
 1974  !! Note that there is no temperature dependence in phytoplankton 
 1975  !! growth rate or any other function. 
 1976  !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid 
 1977  !! NaNs in case of Phy==0. 
 1978  !! 
 1979  !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and nondiat: 
 1980  !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012 
 1981  !! 
 1982  !! AXY (16/07/09) 
 1983  !! temperature for new Eppley style phytoplankton growth 
 1984  loc_T = tsn(ji,jj,jk,jp_tem) 
 1985  fun_T = 1.066**(1.0 * loc_T) 
 1986  !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for 
 1987  !phytoplankton 
 1988  !! growth; remin. unaffected 
 1989  fun_Q10 = xq10**((loc_T  0.0) / 10.0) 
 1990  if (jphy.eq.1) then 
 1991  xvpnT = xvpn * fun_T 
 1992  xvpdT = xvpd * fun_T 
 1993  elseif (jphy.eq.2) then 
 1994  xvpnT = xvpn * fun_Q10 
 1995  xvpdT = xvpd * fun_Q10 
 1996  else 
 1997  xvpnT = xvpn 
 1998  xvpdT = xvpd 
 1999  endif 
 2000  !! 
 2001  !! nondiatoms 
 2002  fchn1 = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
 2003  if (fchn1.GT.rsmall) then 
 2004  fchn = xvpnT / (sqrt(fchn1) + tiny(fchn1)) 
 2005  else 
 2006  fchn = 0. 
 2007  endif 
 2008  fjln = fchn * faln * xpar(ji,jj,jk) !! nondiatom J term 
 2009  fjlim_pn = fjln / xvpnT 
 2010  !! 
 2011  !! diatoms 
 2012  fchd1 = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
 2013  if (fchd1.GT.rsmall) then 
 2014  fchd = xvpdT / (sqrt(fchd1) + tiny(fchd1)) 
 2015  else 
 2016  fchd = 0. 
 2017  endif 
 2018  fjld = fchd * fald * xpar(ji,jj,jk) !! diatom J term 
 2019  fjlim_pd = fjld / xvpdT 
 2020  
 2021  !! 
 2022  !! Phytoplankton nutrient limitation 
 2023  !! 
 2024  !! 
 2025  !! nondiatoms (N, Fe) 
 2026  fnln = zdin / (zdin + xnln) !! nondiatom Qn term 
 2027  ffln = zfer / (zfer + xfln) !! nondiatom Qf term 
 2028  !! 
 2029  !! diatoms (N, Si, Fe) 
 2030  fnld = zdin / (zdin + xnld) !! diatom Qn term 
 2031  fsld = zsil / (zsil + xsld) !! diatom Qs term 
 2032  ffld = zfer / (zfer + xfld) !! diatom Qf term 
 2033  
 2034  !! 
 2035  !! Primary production (nondiatoms) 
 2036  !! (note: still needs multiplying by phytoplankton concentration) 
 2037  !! 
 2038  !! 
 2039  if (jliebig .eq. 0) then 
 2040  !! multiplicative nutrient limitation 
 2041  fpnlim = fnln * ffln 
 2042  elseif (jliebig .eq. 1) then 
 2043  !! Liebig Law (= most limiting) nutrient limitation 
 2044  fpnlim = min(fnln, ffln) 
 2045  endif 
 2046  fprn = fjln * fpnlim 
 2047  
 2048  !! 
 2049  !! Primary production (diatoms) 
 2050  !! (note: still needs multiplying by phytoplankton concentration) 
 2051  !! 
 2052  !! production here is split between nitrogen production and that of 
 2053  !! silicon; depending upon the "intracellular" ratio of Si:N, model 
 2054  !! diatoms will uptake nitrogen/silicon differentially; this borrows 
 2055  !! from the diatom model of Mongin et al. (2006) 
 2056  !! 
 2057  !! 
 2058  if (jliebig .eq. 0) then 
 2059  !! multiplicative nutrient limitation 
 2060  fpdlim = fnld * ffld 
 2061  elseif (jliebig .eq. 1) then 
 2062  !! Liebig Law (= most limiting) nutrient limitation 
 2063  fpdlim = min(fnld, ffld) 
 2064  endif 
 2065  !! 
 2066  if (zphd.GT.rsmall .AND. zpds.GT.rsmall) then 
 2067  !! "intracellular" elemental ratios 
 2068  ! fsin = zpds / (zphd + tiny(zphd)) 
 2069  ! fnsi = zphd / (zpds + tiny(zpds)) 
 2070  fsin = 0.0 
 2071  IF( zphd .GT. rsmall) fsin = zpds / zphd 
 2072  fnsi = 0.0 
 2073  IF( zpds .GT. rsmall) fnsi = zphd / zpds 
 2074  !! AXY (23/02/10): these next variables derive from Mongin et al. (2003) 
 2075  fsin1 = 3.0 * xsin0 !! = 0.6 
 2076  fnsi1 = 1.0 / fsin1 !! = 1.667 
 2077  fnsi2 = 1.0 / xsin0 !! = 5.0 
 2078  !! 
 2079  !! conditionalities based on ratios 
 2080  !! nitrogen (and iron and carbon) 
 2081  if (fsin.le.xsin0) then 
 2082  fprd = 0.0 
 2083  fsld2 = 0.0 
 2084  elseif (fsin.lt.fsin1) then 
 2085  fprd = xuif * ((fsin  xsin0) / (fsin + tiny(fsin))) * (fjld * fpdlim) 
 2086  fsld2 = xuif * ((fsin  xsin0) / (fsin + tiny(fsin))) 
 2087  elseif (fsin.ge.fsin1) then 
 2088  fprd = (fjld * fpdlim) 
 2089  fsld2 = 1.0 
 2090  endif 
 2091  !! 
 2092  !! silicon 
 2093  if (fsin.lt.fnsi1) then 
 2094  fprds = (fjld * fsld) 
 2095  elseif (fsin.lt.fnsi2) then 
 2096  fprds = xuif * ((fnsi  xnsi0) / (fnsi + tiny(fnsi))) * (fjld * fsld) 
 2097  else 
 2098  fprds = 0.0 
 2099  endif 
 2100  else 
 2101  fsin = 0.0 
 2102  fnsi = 0.0 
 2103  fprd = 0.0 
 2104  fsld2 = 0.0 
 2105  fprds = 0.0 
 2106  endif 
 2107  
 2108  !! 
 2109  !! Mixed layer primary production 
 2110  !! this block calculates the amount of primary production that occurs 
 2111  !! within the upper mixed layer; this allows the separate diagnosis 
 2112  !! of "subsurface" primary production; it does assume that short 
 2113  !! term variability in mixed layer depth doesn't mess with things 
 2114  !! though 
 2115  !! 
 2116  !! 
 2117  if (fdep1.le.hmld(ji,jj)) then 
 2118  !! this level is entirely in the mixed layer 
 2119  fq0 = 1.0 
 2120  elseif (fdep.ge.hmld(ji,jj)) then 
 2121  !! this level is entirely below the mixed layer 
 2122  fq0 = 0.0 
 2123  else 
 2124  !! this level straddles the mixed layer 
 2125  fq0 = (hmld(ji,jj)  fdep) / fthk 
 2126  endif 
 2127  !! 
 2128  fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn * zphn * fthk * fq0) 
 2129  fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd * zphd * fthk * fq0) 
 2130  
 2131  !! 
 2132  !! Vertical Integral  
 2133  !! 
 2134  ftot_pn(ji,jj) = ftot_pn(ji,jj) + (zphn * fthk) !! vertical integral nondiatom phytoplankton 
 2135  ftot_pd(ji,jj) = ftot_pd(ji,jj) + (zphd * fthk) !! vertical integral diatom phytoplankton 
 2136  ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi * fthk) !! vertical integral microzooplankton 
 2137  ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme * fthk) !! vertical integral mesozooplankton 
 2138  ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet * fthk) !! vertical integral slow detritus, nitrogen 
 2139  ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc * fthk) !! vertical integral slow detritus, carbon 
 2140  
 2141  !! 
 2142  !! More chlorophyll calculations 
 2143  !! 
 2144  !! 
 2145  !! frn = (xthetam / fthetan) * (fprn / (fthetan * xpar(ji,jj,jk))) 
 2146  !! frd = (xthetam / fthetad) * (fprd / (fthetad * xpar(ji,jj,jk))) 
 2147  frn = (xthetam * fchn * fnln * ffln ) / (fthetan + tiny(fthetan)) 
 2148  !! AXY (12/05/09): there's potentially a problem here; fsld, silicic acid 
 2149  !! limitation, is used in the following line to regulate chlorophyll 
 2150  !! growth in a manner that is inconsistent with its use in the regulation 
 2151  !! of biomass growth; the Mongin term term used in growth is more complex 
 2152  !! than the simple multiplicative function used below 
 2153  !! frd = (xthetam * fchd * fnld * ffld * fsld) / (fthetad + tiny(fthetad)) 
 2154  !! AXY (12/05/09): this replacement line uses the new variable, fsld2, to 
 2155  !! regulate chlorophyll growth 
 2156  frd = (xthetamd * fchd * fnld * ffld * fsld2) / (fthetad + tiny(fthetad)) 
 2157  
 2158  !!====================================================================== 
 2159  !! AXY (07/04/17): possible subroutine block; zooplankton grazing 
 2160  !!====================================================================== 
 2161  
 2162  !! 
 2163  !! Zooplankton Grazing 
 2164  !! this code supplements the base grazing model with one that 
 2165  !! considers the C:N ratio of grazed food and balances this against 
 2166  !! the requirements of zooplankton growth; this model is derived 
 2167  !! from that of Anderson & Pondaven (2003) 
 2168  !! 
 2169  !! the current version of the code assumes a fixed C:N ratio for 
 2170  !! detritus (in contrast to Anderson & Pondaven, 2003), though the 
 2171  !! full equations are retained for future extension 
 2172  !! 
 2173  !! 
 2174  !! 
 2175  !! Microzooplankton first 
 2176  !! 
 2177  !! 
 2178  fmi1 = (xkmi * xkmi) + (xpmipn * zphn * zphn) + (xpmid * zdet * zdet) 
 2179  fmi = xgmi * zzmi / fmi1 
 2180  fgmipn = fmi * xpmipn * zphn * zphn !! grazing on nondiatoms 
 2181  fgmid = fmi * xpmid * zdet * zdet !! grazing on detrital nitrogen 
 2182  # if defined key_roam 
 2183  fgmidc = rsmall !acc 
 2184  IF ( zdet .GT. rsmall ) fgmidc = (zdtc / (zdet + tiny(zdet))) * fgmid !! grazing on detrital carbon 
 2185  # else 
 2186  !! AXY (26/11/08): implicit detrital carbon change 
 2187  fgmidc = xthetad * fgmid !! grazing on detrital carbon 
 2188  # endif 
 2189  !! 
 2190  !! which translates to these incoming N and C fluxes 
 2191  finmi = (1.0  xphi) * (fgmipn + fgmid) 
 2192  ficmi = (1.0  xphi) * ((xthetapn * fgmipn) + fgmidc) 
 2193  !! 
 2194  !! the ideal food C:N ratio for microzooplankton 
 2195  !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
 2196  fstarmi = (xbetan * xthetazmi) / (xbetac * xkc) 
 2197  !! 
 2198  !! process these to determine proportioning of grazed N and C 
 2199  !! (since there is no explicit consideration of respiration, 
 2200  !! only growth and excretion are calculated here) 
 2201  fmith = (ficmi / (finmi + tiny(finmi))) 
 2202  if (fmith.ge.fstarmi) then 
 2203  fmigrow = xbetan * finmi 
 2204  fmiexcr = 0.0 
 2205  else 
 2206  fmigrow = (xbetac * xkc * ficmi) / xthetazmi 
 2207  fmiexcr = ficmi * ((xbetan / (fmith + tiny(fmith)))  ((xbetac * xkc) / xthetazmi)) 
 2208  endif 
 2209  # if defined key_roam 
 2210  fmiresp = (xbetac * ficmi)  (xthetazmi * fmigrow) 
 2211  # endif 
 2212  
 2213  !! 
 2214  !! Mesozooplankton second 
 2215  !! 
 2216  !! 
 2217  fme1 = (xkme * xkme) + (xpmepn * zphn * zphn) + (xpmepd * zphd * zphd) + & 
 2218  (xpmezmi * zzmi * zzmi) + (xpmed * zdet * zdet) 
 2219  fme = xgme * zzme / fme1 
 2220  fgmepn = fme * xpmepn * zphn * zphn !! grazing on nondiatoms 
 2221  fgmepd = fme * xpmepd * zphd * zphd !! grazing on diatoms 
 2222  fgmepds = fsin * fgmepd !! grazing on diatom silicon 
 2223  fgmezmi = fme * xpmezmi * zzmi * zzmi !! grazing on microzooplankton 
 2224  fgmed = fme * xpmed * zdet * zdet !! grazing on detrital nitrogen 
 2225  # if defined key_roam 
 2226  fgmedc = rsmall !acc 
 2227  IF ( zdet .GT. rsmall ) fgmedc = (zdtc / (zdet + tiny(zdet))) * fgmed !! grazing on detrital carbon 
 2228  # else 
 2229  !! AXY (26/11/08): implicit detrital carbon change 
 2230  fgmedc = xthetad * fgmed !! grazing on detrital carbon 
 2231  # endif 
 2232  !! 
 2233  !! which translates to these incoming N and C fluxes 
 2234  finme = (1.0  xphi) * (fgmepn + fgmepd + fgmezmi + fgmed) 
 2235  ficme = (1.0  xphi) * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & 
 2236  (xthetazmi * fgmezmi) + fgmedc) 
 2237  !! 
 2238  !! the ideal food C:N ratio for mesozooplankton 
 2239  !! xbetan = 0.77; xthetaz = 5.625; xbetac = 0.64; xkc = 0.80 
 2240  fstarme = (xbetan * xthetazme) / (xbetac * xkc) 
 2241  !! 
 2242  !! process these to determine proportioning of grazed N and C 
 2243  !! (since there is no explicit consideration of respiration, 
 2244  !! only growth and excretion are calculated here) 
 2245  fmeth = (ficme / (finme + tiny(finme))) 
 2246  if (fmeth.ge.fstarme) then 
 2247  fmegrow = xbetan * finme 
 2248  fmeexcr = 0.0 
 2249  else 
 2250  fmegrow = (xbetac * xkc * ficme) / xthetazme 
 2251  fmeexcr = ficme * ((xbetan / (fmeth + tiny(fmeth)))  ((xbetac * xkc) / xthetazme)) 
 2252  endif 
 2253  # if defined key_roam 
 2254  fmeresp = (xbetac * ficme)  (xthetazme * fmegrow) 
 2255  # endif 
 2256  
 2257  fzmi_i(ji,jj) = fzmi_i(ji,jj) + fthk * ( & 
 2258  fgmipn + fgmid ) 
 2259  fzmi_o(ji,jj) = fzmi_o(ji,jj) + fthk * ( & 
 2260  fmigrow + (xphi * (fgmipn + fgmid)) + fmiexcr + ((1.0  xbetan) * finmi) ) 
 2261  fzme_i(ji,jj) = fzme_i(ji,jj) + fthk * ( & 
 2262  fgmepn + fgmepd + fgmezmi + fgmed ) 
 2263  fzme_o(ji,jj) = fzme_o(ji,jj) + fthk * ( & 
 2264  fmegrow + (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + fmeexcr + ((1.0  xbetan) * finme) ) 
 2265  
 2266  !!====================================================================== 
 2267  !! AXY (07/04/17): possible subroutine block; miscellaneous plankton losses 
 2268  !!====================================================================== 
 2269  
 2270  !! 
 2271  !! Plankton metabolic losses 
 2272  !! Linear loss processes assumed to be metabolic in origin 
 2273  !! 
 2274  !! 
 2275  fdpn2 = xmetapn * zphn 
 2276  fdpd2 = xmetapd * zphd 
 2277  fdpds2 = xmetapd * zpds 
 2278  fdzmi2 = xmetazmi * zzmi 
 2279  fdzme2 = xmetazme * zzme 
 2280  
 2281  !! 
 2282  !! Plankton mortality losses 
 2283  !! EKP (26/02/09): phytoplankton hyperbolic mortality term introduced 
 2284  !! to improve performance in gyres 
 2285  !! 
 2286  !! 
 2287  !! nondiatom phytoplankton 
 2288  if (jmpn.eq.1) fdpn = xmpn * zphn !! linear 
 2289  if (jmpn.eq.2) fdpn = xmpn * zphn * zphn !! quadratic 
 2290  if (jmpn.eq.3) fdpn = xmpn * zphn * & !! hyperbolic 
 2291  (zphn / (xkphn + zphn)) 
 2292  if (jmpn.eq.4) fdpn = xmpn * zphn * & !! sigmoid 
 2293  ((zphn * zphn) / (xkphn + (zphn * zphn))) 
 2294  !! 
 2295  !! diatom phytoplankton 
 2296  if (jmpd.eq.1) fdpd = xmpd * zphd !! linear 
 2297  if (jmpd.eq.2) fdpd = xmpd * zphd * zphd !! quadratic 
 2298  if (jmpd.eq.3) fdpd = xmpd * zphd * & !! hyperbolic 
 2299  (zphd / (xkphd + zphd)) 
 2300  if (jmpd.eq.4) fdpd = xmpd * zphd * & !! sigmoid 
 2301  ((zphd * zphd) / (xkphd + (zphd * zphd))) 
 2302  fdpds = fdpd * fsin 
 2303  !! 
 2304  !! microzooplankton 
 2305  if (jmzmi.eq.1) fdzmi = xmzmi * zzmi !! linear 
 2306  if (jmzmi.eq.2) fdzmi = xmzmi * zzmi * zzmi !! quadratic 
 2307  if (jmzmi.eq.3) fdzmi = xmzmi * zzmi * & !! hyperbolic 
 2308  (zzmi / (xkzmi + zzmi)) 
 2309  if (jmzmi.eq.4) fdzmi = xmzmi * zzmi * & !! sigmoid 
 2310  ((zzmi * zzmi) / (xkzmi + (zzmi * zzmi))) 
 2311  !! 
 2312  !! mesozooplankton 
 2313  if (jmzme.eq.1) fdzme = xmzme * zzme !! linear 
 2314  if (jmzme.eq.2) fdzme = xmzme * zzme * zzme !! quadratic 
 2315  if (jmzme.eq.3) fdzme = xmzme * zzme * & !! hyperbolic 
 2316  (zzme / (xkzme + zzme)) 
 2317  if (jmzme.eq.4) fdzme = xmzme * zzme * & !! sigmoid 
 2318  ((zzme * zzme) / (xkzme + (zzme * zzme))) 
 2319  
 2320  !!====================================================================== 
 2321  !! AXY (07/04/17): possible subroutine block; detritus processes (fuse with later?) 
 2322  !!====================================================================== 
 2323  
 2324  !! 
 2325  !! Detritus remineralisation 
 2326  !! Constant or temperaturedependent 
 2327  !! 
 2328  !! 
 2329  if (jmd.eq.1) then 
 2330  !! temperaturedependent 
 2331  fdd = xmd * fun_T * zdet 
 2332  # if defined key_roam 
 2333  fddc = xmdc * fun_T * zdtc 
 2334  # endif 
 2335  elseif (jmd.eq.2) then 
 2336  !! AXY (16/05/13): add in Q10based parameterisation (def in nmlst) 
 2337  !! temperaturedependent 
 2338  fdd = xmd * fun_Q10 * zdet 
 2339  #if defined key_roam 
 2340  fddc = xmdc * fun_Q10 * zdtc 
 2341  #endif 
 2342  else 
 2343  !! temperatureindependent 
 2344  fdd = xmd * zdet 
 2345  # if defined key_roam 
 2346  fddc = xmdc * zdtc 
 2347  # endif 
 2348  endif 
 2349  !! 
 2350  !! AXY (22/07/09): accelerate detrital remineralisation in the bottom box 
 2351  if ((jk.eq.jmbathy) .and. jsfd.eq.1) then 
 2352  fdd = 1.0 * zdet 
 2353  # if defined key_roam 
 2354  fddc = 1.0 * zdtc 
 2355  # endif 
 2356  endif 
 2357  
 2358  !! 
 2359  !! Detritus addition to benthos 
 2360  !! If activated, slow detritus in the bottom box will enter the 
 2361  !! benthic pool 
 2362  !! 
 2363  !! 
 2364  if ((jk.eq.jmbathy) .and. jorgben.eq.1) then 
 2365  !! this is the BOTTOM OCEAN BOX > into the benthic pool! 
 2366  !! 
 2367  f_sbenin_n(ji,jj) = (zdet * vsed * 86400.) 
 2368  f_sbenin_fe(ji,jj) = (zdet * vsed * 86400. * xrfn) 
 2369  # if defined key_roam 
 2370  f_sbenin_c(ji,jj) = (zdtc * vsed * 86400.) 
 2371  # else 
 2372  f_sbenin_c(ji,jj) = (zdet * vsed * 86400. * xthetad) 
 2373  # endif 
 2374  endif 
 2375  
 2376  !!====================================================================== 
 2377  !! AXY (07/04/17): possible subroutine block; iron chemistry and scavenging 
 2378  !!====================================================================== 
 2379  
 2380  !! 
 2381  !! Iron chemistry and fractionation 
 2382  !! following the Parekh et al. (2004) scheme adopted by the Met. 
 2383  !! Office, Medusa models total iron but considers "free" and 
 2384  !! ligandbound forms for the purposes of scavenging (only "free" 
 2385  !! iron can be scavenged 
 2386  !! 
 2387  !! 
 2388  !! total iron concentration (mmol Fe / m3 > umol Fe / m3) 
 2389  xFeT = zfer * 1.e3 
 2390  !! 
 2391  !! calculate fractionation (based on DiatHadOCC; in turn based on Parekh et al., 2004) 
 2392  xb_coef_tmp = xk_FeL * (xLgT  xFeT)  1.0 
 2393  xb2M4ac = max(((xb_coef_tmp * xb_coef_tmp) + (4.0 * xk_FeL * xLgT)), 0.0) 
 2394  !! 
 2395  !! "free" ligand concentration 
 2396  xLgF = 0.5 * (xb_coef_tmp + (xb2M4ac**0.5)) / xk_FeL 
 2397  !! 
 2398  !! ligandbound iron concentration 
 2399  xFeL = xLgT  xLgF 
 2400  !! 
 2401  !! "free" iron concentration (and convert to mmol Fe / m3) 
 2402  xFeF = (xFeT  xFeL) * 1.e3 
 2403  xFree(ji,jj)= xFeF / (zfer + tiny(zfer)) 
 2404  !! 
 2405  !! scavenging of iron 
 2406  !! AXY (05/04/17): formerly several schemes, now the only appropriate 
 2407  !! one, with all other options returning no scavenging (trc_nam_medusa 
 2408  !! reports on this) 
 2409  !! 
 2410  if (jiron.eq.1) then 
 2411  !! 
 2412  !! Scheme 1: Dutkiewicz et al. (2005) 
 2413  !! This scheme includes a single scavenging term based solely on a 
 2414  !! fixed rate and the availablility of "free" iron 
 2415  !! 
 2416  !! 
 2417  ffescav = xk_sc_Fe * xFeF ! = mmol/m3/d 
 2418  !! 
 2419  !! 
 2420  !! 
 2421  !! Mick's code contains a further (optional) implicit "scavenging" of 
 2422  !! iron that sets an upper bound on "free" iron concentration, and 
 2423  !! essentially caps the concentration of total iron as xFeL + "free" 
 2424  !! iron; since the former is constrained by a fixed total ligand 
 2425  !! concentration (= 1.0 umol/m3), and the latter isn't allowed above 
 2426  !! this upper bound, total iron is constrained to a maximum of ... 
 2427  !! 
 2428  !! xFeL + min(xFeF, 0.3 umol/m3) = 1.0 + 0.3 = 1.3 umol / m3 
 2429  !! 
 2430  !! In Mick's code, the actual value of total iron is reset to this 
 2431  !! sum (i.e. TFe = FeL + Fe'; but Fe' <= 0.3 umol/m3); this isn't 
 2432  !! our favoured approach to tracer updating here (not least because 
 2433  !! of the leapfrog), so here the amount scavenged is augmented by an 
 2434  !! additional amount that serves to drag total iron back towards that 
 2435  !! expected from this limitation on iron concentration ... 
 2436  !! 
 2437  xmaxFeF = min((xFeF * 1.e3), 0.3) ! = umol/m3 
 2438  !! 
 2439  !! Here, the difference between current total Fe and (FeL + Fe') is 
 2440  !! calculated and added to the scavenging flux already calculated 
 2441  !! above ... 
 2442  !! 
 2443  fdeltaFe = (xFeT  (xFeL + xmaxFeF)) * 1.e3 ! = mmol/m3 
 2444  !! 
 2445  !! This assumes that the "excess" iron is dissipated with a time 
 2446  !! scale of 1 day; seems reasonable to me ... (famous last words) 
 2447  !! 
 2448  ffescav = ffescav + fdeltaFe ! = mmol/m3/d 
 2449  !! 
 2450  # if defined key_deep_fe_fix 
 2451  !! AXY (17/01/13) 
 2452  !! stop scavenging for iron concentrations below 0.5 umol / m3 
 2453  !! at depths greater than 1000 m; this aims to end MEDUSA's 
 2454  !! continual loss of iron at depth without impacting things 
 2455  !! at the surface too much; the justification for this is that 
 2456  !! it appears to be what Mick Follows et al. do in their work 
 2457  !! (as evidenced by the iron initial condition they supplied 
 2458  !! me with); to be honest, it looks like Follow et al. do this 
 2459  !! at shallower depths than 1000 m, but I'll stick with this 
 2460  !! for now; I suspect that this seemingly arbitrary approach 
 2461  !! effectively "parameterises" the particlebased scavenging 
 2462  !! rates that other models use (i.e. at depth there are no 
 2463  !! sinking particles, so scavenging stops); it might be fun 
 2464  !! justifying this in a paper though! 
 2465  !! 
 2466  if ((fdep.gt.1000.) .and. (xFeT.lt.0.5)) then 
 2467  ffescav = 0. 
 2468  endif 
 2469  # endif 
 2470  else 
 2471  !! 
 2472  !! No Scheme: you coward! 
 2473  !! This scheme puts its head in the sand and eskews any decision about 
 2474  !! how iron is removed from the ocean; prepare to get deluged in iron 
 2475  !! you fool! 
 2476  !! 
 2477  ffescav = 0. 
 2478  endif 
 2479  
 2480  !! 
 2481  !! Other iron cycle processes 
 2482  !! 
 2483  !! 
 2484  !! aeolian iron deposition 
 2485  if (jk.eq.1) then 
 2486  !! zirondep is in mmolFe / m2 / day 
 2487  !! ffetop is in mmoldissolvedFe / m3 / day 
 2488  ffetop = zirondep(ji,jj) * xfe_sol / fthk 
 2489  else 
 2490  ffetop = 0.0 
 2491  endif 
 2492  !! 
 2493  !! seafloor iron addition 
 2494  !! AXY (10/07/12): amended to only apply sedimentary flux up to ~500 m down 
 2495  !! if (jk.eq.(mbathy(ji,jj)1).AND.jk.lt.i1100) then 
 2496  if ((jk.eq.jmbathy).AND.jk.le.i0500) then 
 2497  !! Moore et al. (2004) cite a coastal California value of 5 umol/m2/d, but adopt a 
 2498  !! global value of 2 umol/m2/d for all areas < 1100 m; here we use this latter value 
 2499  !! but apply it everywhere 
 2500  !! AXY (21/07/09): actually, let's just apply it below 1100 m (levels 137) 
 2501  ffebot = (xfe_sed / fthk) 
 2502  else 
 2503  ffebot = 0.0 
 2504  endif 
 2505  
 2506  !!====================================================================== 
 2507  !! AXY (07/04/17): possible subroutine block; miscellaneous processes (fuse?) 
 2508  !!====================================================================== 
 2509  
 2510  !! 
 2511  !! Miscellaneous 
 2512  !! 
 2513  !! 
 2514  !! diatom frustule dissolution 
 2515  fsdiss = xsdiss * zpds 
 2516  
 2517  !! 
 2518  !! Slow detritus creation 
 2519  !! 
 2520  !! this variable integrates the creation of slow sinking detritus 
 2521  !! to allow the split between fast and slow detritus to be 
 2522  !! diagnosed 
 2523  fslown = fdpn + fdzmi + ((1.0  xfdfrac1) * fdpd) + & 
 2524  ((1.0  xfdfrac2) * fdzme) + ((1.0  xbetan) * (finmi + finme)) 
 2525  !! 
 2526  !! this variable records the slow detrital sinking flux at this 
 2527  !! particular depth; it is used in the output of this flux at 
 2528  !! standard depths in the diagnostic outputs; needs to be 
 2529  !! adjusted from per second to per day because of parameter vsed 
 2530  fslownflux(ji,jj) = zdet * vsed * 86400. 
 2531  # if defined key_roam 
 2532  !! 
 2533  !! and the same for detrital carbon 
 2534  fslowc = (xthetapn * fdpn) + (xthetazmi * fdzmi) + & 
 2535  (xthetapd * (1.0  xfdfrac1) * fdpd) + & 
 2536  (xthetazme * (1.0  xfdfrac2) * fdzme) + & 
 2537  ((1.0  xbetac) * (ficmi + ficme)) 
 2538  !! 
 2539  !! this variable records the slow detrital sinking flux at this 
 2540  !! particular depth; it is used in the output of this flux at 
 2541  !! standard depths in the diagnostic outputs; needs to be 
 2542  !! adjusted from per second to per day because of parameter vsed 
 2543  fslowcflux(ji,jj) = zdtc * vsed * 86400. 
 2544  # endif 
 2545  
 2546  !! 
 2547  !! Nutrient regeneration 
 2548  !! this variable integrates total nitrogen regeneration down the 
 2549  !! watercolumn; its value is stored and output as a 2D diagnostic; 
 2550  !! the corresponding dissolution flux of silicon (from sources 
 2551  !! other than fast detritus) is also integrated; note that, 
 2552  !! confusingly, the linear loss terms from plankton compartments 
 2553  !! are labelled as fdX2 when one might have expected fdX or fdX1 
 2554  !! 
 2555  !! 
 2556  !! nitrogen 
 2557  fregen = (( (xphi * (fgmipn + fgmid)) + & ! messy feeding 
 2558  (xphi * (fgmepn + fgmepd + fgmezmi + fgmed)) + & ! messy feeding 
 2559  fmiexcr + fmeexcr + fdd + & ! excretion + D remin. 
 2560  fdpn2 + fdpd2 + fdzmi2 + fdzme2) * fthk) ! linear mortality 
 2561  !! 
 2562  !! silicon 
 2563  fregensi = (( fsdiss + ((1.0  xfdfrac1) * fdpds) + & ! dissolution + nonlin. mortality 
 2564  ((1.0  xfdfrac3) * fgmepds) + & ! egestion by zooplankton 
 2565  fdpds2) * fthk) ! linear mortality 
 2566  # if defined key_roam 
 2567  !! 
 2568  !! carbon 
 2569  fregenc = (( (xphi * ((xthetapn * fgmipn) + fgmidc)) + & ! messy feeding 
 2570  (xphi * ((xthetapn * fgmepn) + (xthetapd * fgmepd) + & ! messy feeding 
 2571  (xthetazmi * fgmezmi) + fgmedc)) + & ! messy feeding 
 2572  fmiresp + fmeresp + fddc + & ! respiration + D remin. 
 2573  (xthetapn * fdpn2) + (xthetapd * fdpd2) + & ! linear mortality 
 2574  (xthetazmi * fdzmi2) + (xthetazme * fdzme2)) * fthk) ! linear mortality 
 2575  # endif 
 2576  
 2577  
 2578  !!====================================================================== 
 2579  !! AXY (07/04/17): possible subroutine block; fastsinking detritus 
 2580  !!====================================================================== 
 2581  
 2582  !! 
 2583  !! Fastsinking detritus terms 
 2584  !! "local" variables declared so that conservation can be checked; 
 2585  !! the calculated terms are added to the fastsinking flux later on 
 2586  !! only after the flux entering this level has experienced some 
 2587  !! remineralisation 
 2588  !! note: these fluxes need to be scaled by the level thickness 
 2589  !! 
 2590  !! 
 2591  !! nitrogen: diatom and mesozooplankton mortality 
 2592  ftempn = b0 * ((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) 
 2593  !! 
 2594  !! silicon: diatom mortality and grazed diatoms 
 2595  ftempsi = b0 * ((xfdfrac1 * fdpds) + (xfdfrac3 * fgmepds)) 
 2596  !! 
 2597  !! iron: diatom and mesozooplankton mortality 
 2598  ftempfe = b0 * (((xfdfrac1 * fdpd) + (xfdfrac2 * fdzme)) * xrfn) 
 2599  !! 
 2600  !! carbon: diatom and mesozooplankton mortality 
 2601  ftempc = b0 * ((xfdfrac1 * xthetapd * fdpd) + & 
 2602  (xfdfrac2 * xthetazme * fdzme)) 
 2603  !! 
 2604  # if defined key_roam 
 2605  if (jrratio.eq.0) then 
 2606  !! CaCO3: latitudinallybased fraction of total primary production 
 2607  !! absolute latitude of current grid cell 
 2608  flat = abs(gphit(ji,jj)) 
 2609  !! 0.10 at equator; 0.02 at pole 
 2610  fcaco3 = xcaco3a + ((xcaco3b  xcaco3a) * ((90.0  flat) / 90.0)) 
 2611  elseif (jrratio.eq.1) then 
 2612  !! CaCO3: Ridgwell et al. (2007) submodel, version 1 
 2613  !! this uses SURFACE omega calcite to regulate rain ratio 
 2614  if (f_omcal(ji,jj).ge.1.0) then 
 2615  fq1 = (f_omcal(ji,jj)  1.0)**0.81 
 2616  else 
 2617  fq1 = 0. 
 2618  endif 
 2619  fcaco3 = xridg_r0 * fq1 
 2620  elseif (jrratio.eq.2) then 
 2621  !! CaCO3: Ridgwell et al. (2007) submodel, version 2 
 2622  !! this uses FULL 3D omega calcite to regulate rain ratio 
 2623  if (f3_omcal(ji,jj,jk).ge.1.0) then 
 2624  fq1 = (f3_omcal(ji,jj,jk)  1.0)**0.81 
 2625  else 
 2626  fq1 = 0. 
 2627  endif 
 2628  fcaco3 = xridg_r0 * fq1 
 2629  endif 
 2630  # else 
 2631  !! CaCO3: latitudinallybased fraction of total primary production 
 2632  !! absolute latitude of current grid cell 
 2633  flat = abs(gphit(ji,jj)) 
 2634  !! 0.10 at equator; 0.02 at pole 
 2635  fcaco3 = xcaco3a + ((xcaco3b  xcaco3a) * ((90.0  flat) / 90.0)) 
 2636  # endif 
 2637  !! AXY (09/03/09): convert CaCO3 production from function of 
 2638  !! primary production into a function of fastsinking material; 
 2639  !! technically, this is what Dunne et al. (2007) do anyway; they 
 2640  !! convert total primary production estimated from surface 
 2641  !! chlorophyll to an export flux for which they apply conversion 
 2642  !! factors to estimate the various elemental fractions (Si, Ca) 
 2643  ftempca = ftempc * fcaco3 
2077   !! 
2078   !! Phytoplankton light limitation 
2079   !! 
2080   !! 
2081   !! It is assumed xpar is the depthaveraged (vertical layer) PAR 
2082   !! Light limitation (check selfshading) in W/m2 
2083   !! 
2084   !! Note that there is no temperature dependence in phytoplankton 
2085   !! growth rate or any other function. 
2086   !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to avoid 
2087   !! NaNs in case of Phy==0. 
2088   !! 
2089   !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and nondiat: 
2090   !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012 
2091   !! 
2092   !! AXY (16/07/09) 
2093   !! temperature for new Eppley style phytoplankton growth 
2094   loc_T = tsn(ji,jj,jk,jp_tem) 
2095   fun_T = 1.066**(1.0 * loc_T) 
2096   !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for 
2097   !phytoplankton 
2098   !! growth; remin. unaffected 
2099   fun_Q10 = jq10**((loc_T  0.0) / 10.0) 
2100   if (jphy.eq.1) then 
2101   xvpnT = xvpn * fun_T 
2102   xvpdT = xvpd * fun_T 
2103   elseif (jphy.eq.2) then 
2104   xvpnT = xvpn * fun_Q10 
2105   xvpdT = xvpd * fun_Q10 
2106   else 
2107   xvpnT = xvpn 
2108   xvpdT = xvpd 
2109   endif 
2110   !! 
2111   !! nondiatoms 
2112   fchn1 = (xvpnT * xvpnT) + (faln * faln * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
2113   if (fchn1.GT.rsmall) then 
2114   fchn = xvpnT / (sqrt(fchn1) + tiny(fchn1)) 
2115   else 
2116   fchn = 0. 
2117   endif 
2118   fjln = fchn * faln * xpar(ji,jj,jk) !! nondiatom J term 
2119   fjlim_pn = fjln / xvpnT 
2120   !! 
2121   !! diatoms 
2122   fchd1 = (xvpdT * xvpdT) + (fald * fald * xpar(ji,jj,jk) * xpar(ji,jj,jk)) 
2123   if (fchd1.GT.rsmall) then 
2124   fchd = xvpdT / (sqrt(fchd1) + tiny(fchd1)) 
2125   else 
2126   fchd = 0. 
2127   endif 
2128   fjld = fchd * fald * xpar(ji,jj,jk) !! diatom J term 
2129   fjlim_pd = fjld / xvpdT 
2130   
 2657  !! 
 2658  !! This version of MEDUSA offers a choice of three methods for 
 2659  !! handling the remineralisation of fast detritus. All three 
 2660  !! do so in broadly the same way: 
 2661  !! 
 2662  !! 1. Fast detritus is stored as a 2D array [ ffastX ] 
 2663  !! 2. Fast detritus is added levelbylevel [ ftempX ] 
 2664  !! 3. Fast detritus is not remineralised in the top box [ freminX ] 
 2665  !! 4. Remaining fast detritus is remineralised in the bottom [ fsedX ] 
 2666  !! box 
 2667  !! 
 2668  !! The three remineralisation methods are: 
 2669  !! 
 2670  !! 1. Ballast model (i.e. that published in Yool et al., 2011) 
 2671  !! (1b. Ballastsansballast model) 
 2672  !! 2. Martin et al. (1987) 
 2673  !! 3. Henson et al. (2011) 
 2674  !! 
 2675  !! The first of these couples C, N and Fe remineralisation to 
 2676  !! the remineralisation of particulate Si and CaCO3, but the 
 2677  !! latter two treat remineralisation of C, N, Fe, Si and CaCO3 
 2678  !! completely separately. At present a switch within the code 
 2679  !! regulates which submodel is used, but this should be moved 
 2680  !! to the namelist file. 
 2681  !! 
 2682  !! The ballastsansballast submodel is an original development 
 2683  !! feature of MEDUSA in which the ballast submodel's general 
 2684  !! framework and parameterisation is used, but in which there 
 2685  !! is no protection of organic material afforded by ballasting 
 2686  !! minerals. While similar, it is not the same as the Martin 
 2687  !! et al. (1987) submodel. 
 2688  !! 
 2689  !! Since the three submodels behave the same in terms of 
 2690  !! accumulating sinking material and remineralising it all at 
 2691  !! the seafloor, these portions of the code below are common to 
 2692  !! all three. 
 2693  !! 
 2694  
 2695  if (jexport.eq.1) then 
 2696  !!====================================================================== 
 2697  !! BALLAST SUBMODEL 
 2698  !!====================================================================== 
 2699  !! 
 2700  !! 
 2701  !! Fastsinking detritus fluxes, pt. 1: REMINERALISATION 
 2702  !! aside from explicitly modelled, slowsinking detritus, the 
 2703  !! model includes an implicit representation of detrital 
 2704  !! particles that sink too quickly to be modelled with 
 2705  !! explicit state variables; this sinking flux is instead 
 2706  !! instantaneously remineralised down the water column using 
 2707  !! the version of Armstrong et al. (2002)'s ballast model 
 2708  !! used by Dunne et al. (2007); the version of this model 
 2709  !! here considers silicon and calcium carbonate ballast 
 2710  !! minerals; this section of the code redistributes the fast 
 2711  !! sinking material generated locally down the water column; 
 2712  !! this differs from Dunne et al. (2007) in that fast sinking 
 2713  !! material is distributed at *every* level below that it is 
 2714  !! generated, rather than at every level below some fixed 
 2715  !! depth; this scheme is also different in that sinking material 
 2716  !! generated in one level is aggregated with that generated by 
 2717  !! shallower levels; this should make the ballast model more 
 2718  !! selfconsistent (famous last words) 
 2719  !! 
 2720  !! 
 2721  if (jk.eq.1) then 
 2722  !! this is the SURFACE OCEAN BOX (no remineralisation) 
 2723  !! 
 2724  freminc = 0.0 
 2725  freminn = 0.0 
 2726  freminfe = 0.0 
 2727  freminsi = 0.0 
 2728  freminca = 0.0 
 2729  elseif (jk.le.jmbathy) then 
 2730  !! this is an OCEAN BOX (remineralise some material) 
 2731  !! 
 2732  !! set up CCD depth to be used depending on user choice 
 2733  if (jocalccd.eq.0) then 
 2734  !! use default CCD field 
 2735  fccd_dep = ocal_ccd(ji,jj) 
 2736  elseif (jocalccd.eq.1) then 
 2737  !! use calculated CCD field 
 2738  fccd_dep = f2_ccd_cal(ji,jj) 
 2739  endif 
 2740  !! 
 2741  !! === organic carbon === 
 2742  fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) 
 2743  if (iball.eq.1) then 
 2744  fq1 = (fq0 * xmassc) !! how much it weighs (mass) 
 2745  fq2 = (ffastca(ji,jj) * xmassca) !! how much CaCO3 enters this box (mass) 
 2746  fq3 = (ffastsi(ji,jj) * xmasssi) !! how much opal enters this box (mass) 
 2747  fq4 = (fq2 * xprotca) + (fq3 * xprotsi) !! total protected organic C (mass) 
 2748  !! this next term is calculated for C but used for N and Fe as well 
 2749  !! it needs to be protected in case ALL C is protected 
 2750  if (fq4.lt.fq1) then 
 2751  fprotf = (fq4 / (fq1 + tiny(fq1))) !! protected fraction of total organic C (nondim) 
 2752  else 
 2753  fprotf = 1.0 !! all organic C is protected (nondim) 
 2754  endif 
 2755  fq5 = (1.0  fprotf) !! unprotected fraction of total organic C (nondim) 
 2756  fq6 = (fq0 * fq5) !! how much organic C is unprotected (mol) 
 2757  fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected C leaves this box (mol) 
 2758  fq8 = (fq7 + (fq0 * fprotf)) !! how much total C leaves this box (mol) 
 2759  freminc = (fq0  fq8) / fthk !! C remineralisation in this box (mol) 
 2760  ffastc(ji,jj) = fq8 
 2761  else 
 2762  fq1 = fq0 * exp((fthk / xfastc)) !! how much organic C leaves this box (mol) 
 2763  freminc = (fq0  fq1) / fthk !! C remineralisation in this box (mol) 
 2764  ffastc(ji,jj) = fq1 
 2765  endif 
 2766  !! 
 2767  !! === organic nitrogen === 
 2768  fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) 
 2769  if (iball.eq.1) then 
 2770  fq5 = (1.0  fprotf) !! unprotected fraction of total organic N (nondim) 
 2771  fq6 = (fq0 * fq5) !! how much organic N is unprotected (mol) 
 2772  fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected N leaves this box (mol) 
 2773  fq8 = (fq7 + (fq0 * fprotf)) !! how much total N leaves this box (mol) 
 2774  freminn = (fq0  fq8) / fthk !! N remineralisation in this box (mol) 
 2775  ffastn(ji,jj) = fq8 
 2776  else 
 2777  fq1 = fq0 * exp((fthk / xfastc)) !! how much organic N leaves this box (mol) 
 2778  freminn = (fq0  fq1) / fthk !! N remineralisation in this box (mol) 
 2779  ffastn(ji,jj) = fq1 
 2780  endif 
 2781  !! 
 2782  !! === organic iron === 
 2783  fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) 
 2784  if (iball.eq.1) then 
 2785  fq5 = (1.0  fprotf) !! unprotected fraction of total organic Fe (nondim) 
 2786  fq6 = (fq0 * fq5) !! how much organic Fe is unprotected (mol) 
 2787  fq7 = (fq6 * exp((fthk / xfastc))) !! how much unprotected Fe leaves this box (mol) 
 2788  fq8 = (fq7 + (fq0 * fprotf)) !! how much total Fe leaves this box (mol) 
 2789  freminfe = (fq0  fq8) / fthk !! Fe remineralisation in this box (mol) 
 2790  ffastfe(ji,jj) = fq8 
 2791  else 
 2792  fq1 = fq0 * exp((fthk / xfastc)) !! how much total Fe leaves this box (mol) 
 2793  freminfe = (fq0  fq1) / fthk !! Fe remineralisation in this box (mol) 
 2794  ffastfe(ji,jj) = fq1 
 2795  endif 
 2796  !! 
 2797  !! === biogenic silicon === 
 2798  fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) 
 2799  fq1 = fq0 * exp((fthk / xfastsi)) !! how much opal leaves this box (mol) 
 2800  freminsi = (fq0  fq1) / fthk !! Si remineralisation in this box (mol) 
 2801  ffastsi(ji,jj) = fq1 
 2802  !! 
 2803  !! === biogenic calcium carbonate === 
 2804  fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) 
 2805  if (fdep.le.fccd_dep) then 
 2806  !! whole grid cell above CCD 
 2807  fq1 = fq0 !! above lysocline, no Ca dissolves (mol) 
 2808  freminca = 0.0 !! above lysocline, no Ca dissolves (mol) 
 2809  fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) 
 2810  elseif (fdep.ge.fccd_dep) then 
 2811  !! whole grid cell below CCD 
 2812  fq1 = fq0 * exp((fthk / xfastca)) !! how much CaCO3 leaves this box (mol) 
 2813  freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
 2814  else 
 2815  !! partial grid cell below CCD 
 2816  fq2 = fdep1  fccd_dep !! amount of grid cell below CCD (m) 
 2817  fq1 = fq0 * exp((fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) 
 2818  freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
 2819  endif 
 2820  ffastca(ji,jj) = fq1 
 2821  else 
 2822  !! this is BELOW THE LAST OCEAN BOX (do nothing) 
 2823  freminc = 0.0 
 2824  freminn = 0.0 
 2825  freminfe = 0.0 
 2826  freminsi = 0.0 
 2827  freminca = 0.0 
 2828  endif 
 2829  
 2830  elseif (jexport.eq.2.or.jexport.eq.3) then 
 2831  if (jexport.eq.2) then 
 2832  !!====================================================================== 
 2833  !! MARTIN ET AL. (1987) SUBMODEL 
 2834  !!====================================================================== 
 2835  !! 
 2836  !! 
 2837  !! This submodel uses the classic Martin et al. (1987) curve 
 2838  !! to determine the attenuation of fastsinking detritus down 
 2839  !! the water column. All three organic elements, C, N and Fe, 
 2840  !! are handled identically, and their quantities in sinking 
 2841  !! particles attenuate according to a power relationship 
 2842  !! governed by parameter "b". This is assigned a canonical 
 2843  !! value of 0.858. Biogenic opal and calcium carbonate are 
 2844  !! attentuated using the same function as in the ballast 
 2845  !! submodel 
 2846  !! 
 2847  !! 
 2848  fb_val = 0.858 
 2849  elseif (jexport.eq.3) then 
 2850  !!====================================================================== 
 2851  !! HENSON ET AL. (2011) SUBMODEL 
 2852  !!====================================================================== 
 2853  !! 
 2854  !! 
 2855  !! This submodel reconfigures the Martin et al. (1987) curve by 
 2856  !! allowing the "b" value to vary geographically. Its value is 
 2857  !! set, following Henson et al. (2011), as a function of local 
 2858  !! sea surface temperature: 
 2859  !! b = 1.06 + (0.024 * SST) 
 2860  !! This means that remineralisation length scales are longer in 
 2861  !! warm, tropical areas and shorter in cold, polar areas. This 
 2862  !! does seem backtofront (i.e. one would expect GREATER 
 2863  !! remineralisation in warmer waters), but is an outcome of 
 2864  !! analysis of sediment trap data, and it may reflect details 
 2865  !! of ecosystem structure that pertain to particle production 
 2866  !! rather than simply Q10. 
 2867  !! 
 2868  !! 
 2869  fl_sst = tsn(ji,jj,1,jp_tem) 
 2870  fb_val = 1.06 + (0.024 * fl_sst) 
 2871  endif 
 2872  !! 
 2873  if (jk.eq.1) then 
 2874  !! this is the SURFACE OCEAN BOX (no remineralisation) 
 2875  !! 
 2876  freminc = 0.0 
 2877  freminn = 0.0 
 2878  freminfe = 0.0 
 2879  freminsi = 0.0 
 2880  freminca = 0.0 
 2881  elseif (jk.le.jmbathy) then 
 2882  !! this is an OCEAN BOX (remineralise some material) 
 2883  !! 
 2884  !! === organic carbon === 
 2885  fq0 = ffastc(ji,jj) !! how much organic C enters this box (mol) 
 2886  fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic C leaves this box (mol) 
 2887  freminc = (fq0  fq1) / fthk !! C remineralisation in this box (mol) 
 2888  ffastc(ji,jj) = fq1 
 2889  !! 
 2890  !! === organic nitrogen === 
 2891  fq0 = ffastn(ji,jj) !! how much organic N enters this box (mol) 
 2892  fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic N leaves this box (mol) 
 2893  freminn = (fq0  fq1) / fthk !! N remineralisation in this box (mol) 
 2894  ffastn(ji,jj) = fq1 
 2895  !! 
 2896  !! === organic iron === 
 2897  fq0 = ffastfe(ji,jj) !! how much organic Fe enters this box (mol) 
 2898  fq1 = fq0 * ((fdep1/fdep)**fb_val) !! how much organic Fe leaves this box (mol) 
 2899  freminfe = (fq0  fq1) / fthk !! Fe remineralisation in this box (mol) 
 2900  ffastfe(ji,jj) = fq1 
 2901  !! 
 2902  !! === biogenic silicon === 
 2903  fq0 = ffastsi(ji,jj) !! how much opal centers this box (mol) 
 2904  fq1 = fq0 * exp((fthk / xfastsi)) !! how much opal leaves this box (mol) 
 2905  freminsi = (fq0  fq1) / fthk !! Si remineralisation in this box (mol) 
 2906  ffastsi(ji,jj) = fq1 
 2907  !! 
 2908  !! === biogenic calcium carbonate === 
 2909  fq0 = ffastca(ji,jj) !! how much CaCO3 enters this box (mol) 
 2910  if (fdep.le.ocal_ccd(ji,jj)) then 
 2911  !! whole grid cell above CCD 
 2912  fq1 = fq0 !! above lysocline, no Ca dissolves (mol) 
 2913  freminca = 0.0 !! above lysocline, no Ca dissolves (mol) 
 2914  fccd(ji,jj) = real(jk) !! which is the last level above the CCD? (#) 
 2915  elseif (fdep.ge.ocal_ccd(ji,jj)) then 
 2916  !! whole grid cell below CCD 
 2917  fq1 = fq0 * exp((fthk / xfastca)) !! how much CaCO3 leaves this box (mol) 
 2918  freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
 2919  else 
 2920  !! partial grid cell below CCD 
 2921  fq2 = fdep1  ocal_ccd(ji,jj) !! amount of grid cell below CCD (m) 
 2922  fq1 = fq0 * exp((fq2 / xfastca)) !! how much CaCO3 leaves this box (mol) 
 2923  freminca = (fq0  fq1) / fthk !! Ca remineralisation in this box (mol) 
 2924  endif 
 2925  ffastca(ji,jj) = fq1 
 2926  else 
 2927  !! this is BELOW THE LAST OCEAN BOX (do nothing) 
 2928  freminc = 0.0 
 2929  freminn = 0.0 
 2930  freminfe = 0.0 
 2931  freminsi = 0.0 
 2932  freminca = 0.0 
 2933  endif 
 2934  
