- Timestamp:
- 2018-11-13T18:21:16+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
r8427 r10302 30 30 USE daymod 31 31 !! AXY (05/11/13): need these for MEDUSA to input/output benthic reservoirs 32 USE par_medusa 32 33 USE sms_medusa 33 34 USE trcsms_medusa … … 43 44 USE sbc_oce, ONLY: lk_oasis 44 45 USE oce, ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl !! Coupling variable 46 USE trcstat 47 USE obs_const, ONLY: obfillflt ! Observation operator fill value 45 48 46 49 IMPLICIT NONE … … 52 55 PUBLIC trc_rst_cal 53 56 PUBLIC trc_rst_stat 54 PUBLIC trc_rst_dia_stat 55 PUBLIC trc_rst_tra_stat 57 #if defined key_medusa && defined key_roam 58 PUBLIC trc_rst_conserve 59 #endif 56 60 57 61 !! * Substitutions … … 277 281 !! as proxy of org matter from the ocean 278 282 !! -- needed for the coupling with atm 283 !! 07-12-2017 -- To make things cleaner, we want to store an 284 !! unscaled Chl field in the restart and only 285 !! scale it when reading it in. 286 279 287 IF( iom_varid( numrtr, 'N_CHL_srf', ldstop = .FALSE. ) > 0 ) THEN 280 IF(lwp) WRITE(numout,*) 'Chl surf concentration - reading in ...'288 IF(lwp) WRITE(numout,*) 'Chl cpl concentration - reading in ... - scale by ', scl_chl 281 289 CALL iom_get( numrtr, jpdom_autoglo, 'N_CHL_srf', zn_chl_srf(:,:) ) 282 290 ELSE 283 IF(lwp) WRITE(numout,*) ' Chl surf concentration - setting to zero ...'284 zn_chl_srf(:,:) = (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6291 IF(lwp) WRITE(numout,*) 'set Chl coupled concentration - scaled by ', scl_chl 292 zn_chl_srf(:,:) = MAX( 0.0, (trn(:,:,1,jpchn) + trn(:,:,1,jpchd)) * 1.E-6 ) 285 293 ENDIF 286 294 IF (lk_oasis) THEN 287 chloro_out_cpl(:,:) = zn_chl_srf(:,:) !! Coupling variable295 chloro_out_cpl(:,:) = zn_chl_srf(:,:) * scl_chl !! Coupling variable 288 296 END IF 289 297 !! … … 297 305 call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf') 298 306 call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux') 299 call trc_rst_dia_stat(zn_chl_srf(:,:), 'CHL surf') 307 IF (lk_oasis) THEN 308 call trc_rst_dia_stat(chloro_out_cpl(:,:), 'CHL cpl') 309 END IF 300 310 !! 301 311 !! JPALM 14-06-2016 -- add Carbonate chenistry variables through the restart … … 329 339 IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 330 340 ENDIF 341 ! 342 IF ( ln_foam_medusa ) THEN 343 !! 2D fields of pCO2 and fCO2 for observation operator on first timestep 344 IF( iom_varid( numrtr, 'PCO2W', ldstop = .FALSE. ) > 0 ) THEN 345 IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 present - reading in ...' 346 CALL iom_get( numrtr, jpdom_autoglo, 'PCO2W', f2_pco2w(:,:) ) 347 CALL iom_get( numrtr, jpdom_autoglo, 'FCO2W', f2_fco2w(:,:) ) 348 ELSE 349 IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 absent - setting to fill ...' 350 f2_pco2w(:,:) = obfillflt * tmask(:,:,1) 351 f2_fco2w(:,:) = obfillflt * tmask(:,:,1) 352 ENDIF 353 ENDIF 331 354 # endif 332 355 IF ( ln_foam_medusa ) THEN 356 !! Fields for ocean colour assimilation on first timestep 357 IF( iom_varid( numrtr, 'pgrow_avg', ldstop = .FALSE. ) > 0 ) THEN 358 IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...' 359 CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg', pgrow_avg(:,:) ) 360 CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg', ploss_avg(:,:) ) 361 CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg', phyt_avg(:,:) ) 362 CALL iom_get( numrtr, jpdom_autoglo, 'mld_max', mld_max(:,:) ) 363 ELSE 364 IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg absent - setting to zero ...' 365 pgrow_avg(:,:) = 0.0 366 ploss_avg(:,:) = 0.0 367 phyt_avg(:,:) = 0.0 368 mld_max(:,:) = 0.0 369 ENDIF 370 ENDIF 333 371 334 372 #endif … … 369 407 !! 370 408 IF(lwp) WRITE(numout,*) ' CFC averaged properties absent - setting to zero ...' 371 qint_cfc(:,:,j n) = 0.0 !! CHN409 qint_cfc(:,:,jl) = 0.0 !! CHN 372 410 ENDIF 373 411 !! … … 457 495 CALL iom_rstput( kt, nitrst, numrtw, 'B_CO2_flx', zb_co2_flx(:,:) ) 458 496 CALL iom_rstput( kt, nitrst, numrtw, 'N_CO2_flx', zn_co2_flx(:,:) ) 497 !! JPALM 07-12-2017 -- To make things cleaner, we want to store an 498 !! unscaled Chl field in the restart and only 499 !! scale it when reading it in. 459 500 CALL iom_rstput( kt, nitrst, numrtw, 'N_CHL_srf', zn_chl_srf(:,:) ) 460 501 !! … … 468 509 call trc_rst_dia_stat(zn_dms_srf(:,:), 'DMS surf') 469 510 call trc_rst_dia_stat(zn_co2_flx(:,:), 'CO2 flux') 470 call trc_rst_dia_stat(zn_chl_srf(:,:), ' CHL surf')511 call trc_rst_dia_stat(zn_chl_srf(:,:), 'unscaled CHL cpl') 471 512 !! 472 513 IF(lwp) WRITE(numout,*) ' MEDUSA averaged prop. for dust and iron dep.' … … 498 539 call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG') 499 540 !! 541 IF ( ln_foam_medusa ) THEN 542 !! Fields for observation operator on first timestep 543 IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...' 544 CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:) ) 545 CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:) ) 546 ENDIF 500 547 # endif 548 IF ( ln_foam_medusa ) THEN 549 !! Fields for assimilation on first timestep 550 IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...' 551 CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) ) 552 CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) ) 553 CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg', phyt_avg(:,:) ) 554 CALL iom_rstput( kt, nitrst, numrtw, 'mld_max', mld_max(:,:) ) 555 ENDIF 501 556 !! 502 557 #endif … … 531 586 IF( kt == nitrst ) THEN 532 587 CALL trc_rst_stat ! statistics 588 #if defined key_medusa && defined key_roam 589 CALL trc_rst_conserve ! conservation check 590 #endif 533 591 CALL iom_close( numrtw ) ! close the restart file (only at last time step) 534 592 #if ! defined key_trdmxl_trc … … 697 755 698 756 699 SUBROUTINE trc_rst_tra_stat 700 !!---------------------------------------------------------------------- 701 !! *** trc_rst_tra_stat *** 702 !! 703 !! ** purpose : Compute tracers statistics - check where crazy values appears 704 !!---------------------------------------------------------------------- 705 INTEGER :: jk, jn 706 REAL(wp) :: ztraf, zmin, zmax, zmean, zdrift, areasf 707 REAL(wp), DIMENSION(jpi,jpj) :: zvol 708 !!---------------------------------------------------------------------- 709 757 # if defined key_medusa && defined key_roam 758 SUBROUTINE trc_rst_conserve 759 !!---------------------------------------------------------------------- 760 !! *** trc_rst_conserve *** 761 !! 762 !! ** purpose : Compute tracers conservation statistics 763 !! 764 !! AXY (17/11/2017) 765 !! This routine calculates the "now" inventories of the elemental 766 !! cycles of MEDUSA and compares them to those calculate when the 767 !! model was initialised / restarted; the cycles calculated are: 768 !! nitrogen, silicon, iron, carbon, alkalinity and oxygen 769 !!---------------------------------------------------------------------- 770 INTEGER :: ji, jj, jk, jn 771 REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio 772 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol 773 REAL(wp), DIMENSION(jpi,jpj) :: z2d, zarea 774 REAL(wp), DIMENSION(6) :: loc_cycletot3, loc_cycletot2 775 !!---------------------------------------------------------------------- 776 ! 710 777 IF( lwp ) THEN 778 WRITE(numout,*) 779 WRITE(numout,*) ' ----TRACER CONSERVATION---- ' 780 WRITE(numout,*) 781 ENDIF 782 ! 783 ! ocean volume 784 DO jk = 1, jpk 785 zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk) 786 END DO 787 ! 788 ! ocean area (for sediments) 789 zarea(:,:) = e1e2t(:,:) * tmask(:,:,1) 790 ! 791 !---------------------------------------------------------------------- 792 ! nitrogen 793 z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 794 trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 795 z2d(:,:) = zn_sed_n(:,:) 796 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 797 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 798 ! total tracer, and delta 799 zinvt = zsum3d + zsum2d 800 zdelta = zinvt - cycletot(1) 801 zratio = 1.0e2 * zdelta / cycletot(1) 802 ! 803 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt, & 804 cycletot(1), zdelta, zratio 805 IF ( lwp ) WRITE(numout,*) 806 ! 807 !---------------------------------------------------------------------- 808 ! silicon 809 z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil) 810 z2d(:,:) = zn_sed_si(:,:) 811 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 812 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 813 ! total tracer, and delta 814 zinvt = zsum3d + zsum2d 815 zdelta = zinvt - cycletot(2) 816 zratio = 1.0e2 * zdelta / cycletot(2) 817 ! 818 IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt, & 819 cycletot(2), zdelta, zratio 820 IF ( lwp ) WRITE(numout,*) 821 ! 822 !---------------------------------------------------------------------- 823 ! iron 824 z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 825 trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer) 826 z2d(:,:) = zn_sed_fe(:,:) 827 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 828 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 829 ! total tracer, and delta 830 zinvt = zsum3d + zsum2d 831 zdelta = zinvt - cycletot(3) 832 zratio = 1.0e2 * zdelta / cycletot(3) 833 ! 834 IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt, & 835 cycletot(3), zdelta, zratio 836 IF ( lwp ) WRITE(numout,*) 837 ! 838 !---------------------------------------------------------------------- 839 ! carbon 840 z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn) + (trn(:,:,:,jpphd) * xthetapd) + & 841 (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + & 842 trn(:,:,:,jpdtc) + trn(:,:,:,jpdic) 843 z2d(:,:) = zn_sed_c(:,:) + zn_sed_ca(:,:) 844 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 845 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 846 ! total tracer, and delta 847 zinvt = zsum3d + zsum2d 848 zdelta = zinvt - cycletot(4) 849 zratio = 1.0e2 * zdelta / cycletot(4) 850 ! 851 IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt, & 852 cycletot(4), zdelta, zratio 853 IF ( lwp ) WRITE(numout,*) 854 ! 855 !---------------------------------------------------------------------- 856 ! alkalinity 857 z3d(:,:,:) = trn(:,:,:,jpalk) 858 z2d(:,:) = zn_sed_ca(:,:) * 2.0 859 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 860 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 861 ! total tracer, and delta 862 zinvt = zsum3d + zsum2d 863 zdelta = zinvt - cycletot(5) 864 zratio = 1.0e2 * zdelta / cycletot(5) 865 ! 866 IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, & 867 cycletot(5), zdelta, zratio 868 IF ( lwp ) WRITE(numout,*) 869 ! 870 !---------------------------------------------------------------------- 871 ! oxygen 872 z3d(:,:,:) = trn(:,:,:,jpoxy) 873 z2d(:,:) = 0.0 874 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 875 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 876 ! total tracer, and delta 877 zinvt = zsum3d + zsum2d 878 zdelta = zinvt - cycletot(6) 879 zratio = 1.0e2 * zdelta / cycletot(6) 880 ! 881 IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt, & 882 cycletot(6), zdelta, zratio 883 ! 884 !---------------------------------------------------------------------- 885 ! Check 886 zsum3d = glob_sum( zvol(:,:,:) ) 887 zsum2d = glob_sum( zarea(:,:) ) 888 IF ( lwp ) THEN 711 889 WRITE(numout,*) 712 WRITE(numout,*) ' ----SURFACE TRA STAT---- ' 890 WRITE(numout,*) ' check : cvol : ', zsum3d 891 WRITE(numout,*) ' check : carea : ', zsum2d 713 892 WRITE(numout,*) 714 893 ENDIF 715 894 ! 716 zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1) 717 areasf = glob_sum(zvol(:,:)) 718 DO jn = 1, jptra 719 ztraf = glob_sum( tra(:,:,1,jn) * zvol(:,:) ) 720 zmin = MINVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) ) 721 zmax = MAXVAL( tra(:,:,1,jn), mask= ((tmask(:,:,1).NE.0.)) ) 722 IF( lk_mpp ) THEN 723 CALL mpp_min( zmin ) ! min over the global domain 724 CALL mpp_max( zmax ) ! max over the global domain 725 END IF 726 zmean = ztraf / areasf 727 IF(lwp) WRITE(numout,9001) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax 728 END DO 729 IF(lwp) WRITE(numout,*) 730 9001 FORMAT(' tracer nb :',i2,' name :',a10,' mean :',e18.10,' min :',e18.10, & 731 & ' max :',e18.10) 732 ! 733 END SUBROUTINE trc_rst_tra_stat 734 735 736 737 SUBROUTINE trc_rst_dia_stat( dgtr, names) 738 !!---------------------------------------------------------------------- 739 !! *** trc_rst_dia_stat *** 740 !! 741 !! ** purpose : Compute tracers statistics 742 !!---------------------------------------------------------------------- 743 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: dgtr ! 2D diag var 744 CHARACTER(len=*) , INTENT(in) :: names ! 2D diag name 745 !!--------------------------------------------------------------------- 746 INTEGER :: jk, jn 747 CHARACTER (LEN=18) :: text_zmean 748 REAL(wp) :: ztraf, zmin, zmax, zmean, areasf 749 REAL(wp), DIMENSION(jpi,jpj) :: zvol 750 !!---------------------------------------------------------------------- 751 752 IF( lwp ) WRITE(numout,*) 'STAT- ', names 753 754 ! fse3t_a will be undefined at the start of a run, but this routine 755 ! may be called at any stage! Hence we MUST make sure it is 756 ! initialised to zero when allocated to enable us to test for 757 ! zero content here and avoid potentially dangerous and non-portable 758 ! operations (e.g. divide by zero, global sums of junk values etc.) 759 zvol(:,:) = e1e2t(:,:) * fse3t_a(:,:,1) * tmask(:,:,1) 760 ztraf = glob_sum( dgtr(:,:) * zvol(:,:) ) 761 !! areasf = glob_sum(e1e2t(:,:) * tmask(:,:,1) ) 762 areasf = glob_sum(zvol(:,:)) 763 zmin = MINVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) ) 764 zmax = MAXVAL( dgtr(:,:), mask= ((tmask(:,:,1).NE.0.)) ) 765 IF( lk_mpp ) THEN 766 CALL mpp_min( zmin ) ! min over the global domain 767 CALL mpp_max( zmax ) ! max over the global domain 768 END IF 769 770 text_zmean = "N/A" 771 ! Avoid divide by zero. areasf must be positive. 772 IF (areasf > 0.0) THEN 773 zmean = ztraf / areasf 774 WRITE(text_zmean,'(e18.10)') zmean 775 ENDIF 776 777 IF(lwp) WRITE(numout,9002) TRIM( names ), text_zmean, zmin, zmax 778 779 9002 FORMAT(' tracer name :',A,' mean :',A,' min :',e18.10, & 780 & ' max :',e18.10 ) 781 ! 782 END SUBROUTINE trc_rst_dia_stat 895 9010 FORMAT(' element:',a10, & 896 ' 3d sum:',e18.10,' 2d sum:',e18.10, & 897 ' total:',e18.10,' initial:',e18.10, & 898 ' delta:',e18.10,' %:',e18.10) 899 ! 900 END SUBROUTINE trc_rst_conserve 901 # endif 783 902 784 903
Note: See TracChangeset
for help on using the changeset viewer.
