- Timestamp:
- 2017-12-06T17:10:46+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5518_GO6_conserv_Check/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
r8427 r8926 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 … … 52 53 PUBLIC trc_rst_cal 53 54 PUBLIC trc_rst_stat 55 #if defined key_medusa && defined key_roam 56 PUBLIC trc_rst_conserve 57 #endif 54 58 PUBLIC trc_rst_dia_stat 55 59 PUBLIC trc_rst_tra_stat … … 531 535 IF( kt == nitrst ) THEN 532 536 CALL trc_rst_stat ! statistics 537 #if defined key_medusa && defined key_roam 538 CALL trc_rst_conserve ! conservation check 539 #endif 533 540 CALL iom_close( numrtw ) ! close the restart file (only at last time step) 534 541 #if ! defined key_trdmxl_trc … … 697 704 698 705 706 #if defined key_medusa && defined key_roam 707 SUBROUTINE trc_rst_conserve 708 !!---------------------------------------------------------------------- 709 !! *** trc_rst_conserve *** 710 !! 711 !! ** purpose : Compute tracers conservation statistics 712 !! 713 !! AXY (17/11/2017) 714 !! This routine calculates the "now" inventories of the elemental 715 !! cycles of MEDUSA and compares them to those calculate when the 716 !! model was initialised / restarted; the cycles calculated are: 717 !! nitrogen, silicon, iron, carbon, alkalinity and oxygen 718 !!---------------------------------------------------------------------- 719 INTEGER :: ji, jj, jk, jn 720 REAL(wp) :: zsum3d, zsum2d, zinvt, zdelta, zratio, loc_vol, loc_are 721 REAL(wp) :: zq1, zq2, loc_vol, loc_area 722 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d, zvol 723 REAL(wp), DIMENSION(jpi,jpj) :: z2d, zarea 724 REAL(wp), DIMENSION(6) :: loc_cycletot3, loc_cycletot2 725 !!---------------------------------------------------------------------- 726 ! 727 IF( lwp ) THEN 728 WRITE(numout,*) 729 WRITE(numout,*) ' ----TRACER CONSERVATION---- ' 730 WRITE(numout,*) 731 ENDIF 732 ! 733 ! ocean volume 734 DO jk = 1, jpk 735 zvol(:,:,jk) = e1e2t(:,:) * fse3t_a(:,:,jk) * tmask(:,:,jk) 736 END DO 737 ! 738 ! ocean area (for sediments) 739 zarea(:,:) = e1e2t(:,:) * tmask(:,:,1) 740 ! 741 !---------------------------------------------------------------------- 742 ! nitrogen 743 z3d(:,:,:) = trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 744 trn(:,:,:,jpzme) + trn(:,:,:,jpdet) + trn(:,:,:,jpdin) 745 z2d(:,:) = zn_sed_n(:,:) 746 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 747 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 748 ! total tracer, and delta 749 zinvt = zsum3d + zsum2d 750 zdelta = zinvt - cycletot(1) 751 zratio = 1.0e2 * zdelta / cycletot(1) 752 ! 753 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zsum3d, zsum2d, zinvt, & 754 cycletot(1), zdelta, zratio 755 IF ( lwp ) WRITE(numout,*) 756 ! 757 !---------------------------------------------------------------------- 758 ! silicon 759 z3d(:,:,:) = trn(:,:,:,jppds) + trn(:,:,:,jpsil) 760 z2d(:,:) = zn_sed_si(:,:) 761 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 762 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 763 ! total tracer, and delta 764 zinvt = zsum3d + zsum2d 765 zdelta = zinvt - cycletot(2) 766 zratio = 1.0e2 * zdelta / cycletot(2) 767 ! 768 IF ( lwp ) WRITE(numout,9010) 'silicon', zsum3d, zsum2d, zinvt, & 769 cycletot(2), zdelta, zratio 770 IF ( lwp ) WRITE(numout,*) 771 ! 772 !---------------------------------------------------------------------- 773 ! iron 774 z3d(:,:,:) = ((trn(:,:,:,jpphn) + trn(:,:,:,jpphd) + trn(:,:,:,jpzmi) + & 775 trn(:,:,:,jpzme) + trn(:,:,:,jpdet)) * xrfn) + trn(:,:,:,jpfer) 776 z2d(:,:) = zn_sed_fe(:,:) 777 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 778 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 779 ! total tracer, and delta 780 zinvt = zsum3d + zsum2d 781 zdelta = zinvt - cycletot(3) 782 zratio = 1.0e2 * zdelta / cycletot(3) 783 ! 784 IF ( lwp ) WRITE(numout,9010) 'iron', zsum3d, zsum2d, zinvt, & 785 cycletot(3), zdelta, zratio 786 IF ( lwp ) WRITE(numout,*) 787 ! 788 !---------------------------------------------------------------------- 789 ! carbon 790 z3d(:,:,:) = (trn(:,:,:,jpphn) * xthetapn) + (trn(:,:,:,jpphd) * xthetapd) + & 791 (trn(:,:,:,jpzmi) * xthetazmi) + (trn(:,:,:,jpzme) * xthetazme) + & 792 trn(:,:,:,jpdtc) + trn(:,:,:,jpdic) 793 z2d(:,:) = zn_sed_c(:,:) + zn_sed_ca(:,:) 794 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 795 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 796 ! total tracer, and delta 797 zinvt = zsum3d + zsum2d 798 zdelta = zinvt - cycletot(4) 799 zratio = 1.0e2 * zdelta / cycletot(4) 800 ! 801 IF ( lwp ) WRITE(numout,9010) 'carbon', zsum3d, zsum2d, zinvt, & 802 cycletot(4), zdelta, zratio 803 IF ( lwp ) WRITE(numout,*) 804 ! 805 !---------------------------------------------------------------------- 806 ! alkalinity 807 z3d(:,:,:) = trn(:,:,:,jpalk) 808 z2d(:,:) = zn_sed_ca(:,:) * 2.0 809 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 810 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 811 ! total tracer, and delta 812 zinvt = zsum3d + zsum2d 813 zdelta = zinvt - cycletot(5) 814 zratio = 1.0e2 * zdelta / cycletot(5) 815 ! 816 IF ( lwp ) WRITE(numout,9010) 'alkalinity', zsum3d, zsum2d, zinvt, & 817 cycletot(5), zdelta, zratio 818 IF ( lwp ) WRITE(numout,*) 819 ! 820 !---------------------------------------------------------------------- 821 ! oxygen 822 z3d(:,:,:) = trn(:,:,:,jpoxy) 823 z2d(:,:) = 0.0 824 zsum3d = glob_sum( z3d(:,:,:) * zvol(:,:,:) ) 825 zsum2d = glob_sum( z2d(:,:) * zarea(:,:) ) 826 ! total tracer, and delta 827 zinvt = zsum3d + zsum2d 828 zdelta = zinvt - cycletot(6) 829 zratio = 1.0e2 * zdelta / cycletot(6) 830 ! 831 IF ( lwp ) WRITE(numout,9010) 'oxygen', zsum3d, zsum2d, zinvt, & 832 cycletot(6), zdelta, zratio 833 ! 834 !---------------------------------------------------------------------- 835 ! Check 836 zsum3d = glob_sum( zvol(:,:,:) ) 837 zsum2d = glob_sum( zarea(:,:) ) 838 IF ( lwp ) WRITE(numout,*) 839 IF ( lwp ) WRITE(numout,*) ' check : cvol : ', zsum3d 840 IF ( lwp ) WRITE(numout,*) ' check : carea : ', zsum2d 841 IF ( lwp ) WRITE(numout,*) 842 IF ( lwp ) CALL flush(numout) 843 ! 844 !---------------------------------------------------------------------- 845 ! Excluding Halo version: 846 ! 847 ! This second check dodges around the halo points in the grid 848 ! to check that glob_sum is doing what it's supposed to be 849 ! doing; note that the loop ordering here is the "correct" way 850 ! (according to Dr. Google) 851 ! 852 IF ( lwp ) WRITE(numout,*) 853 IF ( lwp ) WRITE(numout,*) ' Elemental cycle totals (check): ' 854 ! ocean cells 855 loc_cycletot3(:) = 0._wp 856 loc_vol = 0._wp 857 DO jk = 1, jpk 858 DO jj = nldj,nlej 859 DO ji = nldi,nlei 860 loc_vol = loc_vol + cvol(ji,jj,jk) 861 ! nitrogen 862 zq1 = trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 863 trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet) + trn(ji,jj,jk,jpdin) 864 loc_cycletot3(1) = loc_cycletot3(1) + ( zq1 * cvol(ji,jj,jk) ) 865 ! silicon 866 zq1 = trn(ji,jj,jk,jppds) + trn(ji,jj,jk,jpsil) 867 loc_cycletot3(2) = loc_cycletot3(2) + ( zq1 * cvol(ji,jj,jk) ) 868 ! iron 869 zq1 = ((trn(ji,jj,jk,jpphn) + trn(ji,jj,jk,jpphd) + trn(ji,jj,jk,jpzmi) + & 870 trn(ji,jj,jk,jpzme) + trn(ji,jj,jk,jpdet)) * xrfn) + trn(ji,jj,jk,jpfer) 871 loc_cycletot3(3) = loc_cycletot3(3) + ( zq1 * cvol(ji,jj,jk) ) 872 ! carbon 873 zq1 = (trn(ji,jj,jk,jpphn) * xthetapn) + (trn(ji,jj,jk,jpphd) * xthetapd) + & 874 (trn(ji,jj,jk,jpzmi) * xthetazmi) + (trn(ji,jj,jk,jpzme) * xthetazme) + & 875 trn(ji,jj,jk,jpdtc) + trn(ji,jj,jk,jpdic) 876 loc_cycletot3(4) = loc_cycletot3(4) + ( zq1 * cvol(ji,jj,jk) ) 877 ! alkalinity 878 zq1 = trn(ji,jj,jk,jpalk) 879 loc_cycletot3(5) = loc_cycletot3(5) + ( zq1 * cvol(ji,jj,jk) ) 880 ! oxygen 881 zq1 = trn(ji,jj,jk,jpoxy) 882 loc_cycletot3(6) = loc_cycletot3(6) + ( zq1 * cvol(ji,jj,jk) ) 883 ENDDO 884 ENDDO 885 ENDDO 886 ! 887 ! sediment cells 888 loc_cycletot2(:) = 0._wp 889 loc_area = 0._wp 890 jk = 1 891 DO jj = nldj,nlej 892 DO ji = nldi,nlei 893 loc_area = loc_area + (e1e2t(ji,jj) * tmask(ji,jj,jk)) 894 ! nitrogen 895 zq1 = zn_sed_n(ji,jj) 896 loc_cycletot2(1) = loc_cycletot2(1) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 897 ! silicon 898 zq1 = zn_sed_si(ji,jj) 899 loc_cycletot2(2) = loc_cycletot2(2) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 900 ! iron 901 zq1 = zn_sed_fe(ji,jj) 902 loc_cycletot2(3) = loc_cycletot2(3) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 903 ! carbon 904 zq1 = zn_sed_c(ji,jj) + zn_sed_ca(ji,jj) 905 loc_cycletot2(4) = loc_cycletot2(4) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 906 ! alkalinity 907 zq1 = zn_sed_ca(ji,jj) * 2._wp 908 loc_cycletot2(5) = loc_cycletot2(5) + (zq1 * e1e2t(ji,jj) * tmask(ji,jj,jk)) 909 ! skip oxygen 910 ENDDO 911 ENDDO 912 IF ( lwp ) WRITE(numout,*) '-- New Check --' 913 ! 914 !---------------------------------------------------------------------- 915 ! nitrogen 916 ! total tracer, and delta 917 zq1 = loc_cycletot3(1) 918 zq2 = loc_cycletot2(1) 919 IF( lk_mpp ) CALL mpp_sum( zq1 ) 920 IF( lk_mpp ) CALL mpp_sum( zq2 ) 921 zinvt = zq1 + zq2 922 zdelta = zinvt - cycletot2(1) 923 zratio = 1.0e2 * zdelta / cycletot2(1) 924 ! 925 IF ( lwp ) WRITE(numout,9010) 'nitrogen', zq1, zq2, zinvt, & 926 cycletot2(1), zdelta, zratio 927 IF ( lwp ) WRITE(numout,*) 928 ! 929 !---------------------------------------------------------------------- 930 ! silicon 931 zq1 = loc_cycletot3(2) 932 zq2 = loc_cycletot2(2) 933 IF( lk_mpp ) CALL mpp_sum( zq1 ) 934 IF( lk_mpp ) CALL mpp_sum( zq2 ) 935 zinvt = zq1 + zq2 936 zdelta = zinvt - cycletot2(2) 937 zratio = 1.0e2 * zdelta / cycletot2(2) 938 ! 939 IF ( lwp ) WRITE(numout,9010) 'silicon', zq1, zq2, zinvt, & 940 cycletot2(2), zdelta, zratio 941 IF ( lwp ) WRITE(numout,*) 942 ! 943 !---------------------------------------------------------------------- 944 ! iron 945 zq1 = loc_cycletot3(3) 946 zq2 = loc_cycletot2(3) 947 IF( lk_mpp ) CALL mpp_sum( zq1 ) 948 IF( lk_mpp ) CALL mpp_sum( zq2 ) 949 zinvt = zq1 + zq2 950 zdelta = zinvt - cycletot2(3) 951 zratio = 1.0e2 * zdelta / cycletot2(3) 952 ! 953 IF ( lwp ) WRITE(numout,9010) 'iron', zq1, zq2, zinvt, & 954 cycletot2(3), zdelta, zratio 955 IF ( lwp ) WRITE(numout,*) 956 ! 957 !---------------------------------------------------------------------- 958 ! carbon 959 zq1 = loc_cycletot3(4) 960 zq2 = loc_cycletot2(4) 961 IF( lk_mpp ) CALL mpp_sum( zq1 ) 962 IF( lk_mpp ) CALL mpp_sum( zq2 ) 963 zinvt = zq1 + zq2 964 zdelta = zinvt - cycletot2(4) 965 zratio = 1.0e2 * zdelta / cycletot2(4) 966 ! 967 IF ( lwp ) WRITE(numout,9010) 'carbon', zq1, zq2, zinvt, & 968 cycletot2(4), zdelta, zratio 969 IF ( lwp ) WRITE(numout,*) 970 ! 971 !---------------------------------------------------------------------- 972 ! alkalinity 973 zq1 = loc_cycletot3(5) 974 zq2 = loc_cycletot2(5) 975 IF( lk_mpp ) CALL mpp_sum( zq1 ) 976 IF( lk_mpp ) CALL mpp_sum( zq2 ) 977 zinvt = zq1 + zq2 978 zdelta = zinvt - cycletot2(5) 979 zratio = 1.0e2 * zdelta / cycletot2(5) 980 ! 981 IF ( lwp ) WRITE(numout,9010) 'alkalinity', zq1, zq2, zinvt, & 982 cycletot2(5), zdelta, zratio 983 IF ( lwp ) WRITE(numout,*) 984 ! 985 !---------------------------------------------------------------------- 986 ! oxygen 987 zq1 = loc_cycletot3(6) 988 zq2 = loc_cycletot2(6) 989 IF( lk_mpp ) CALL mpp_sum( zq1 ) 990 IF( lk_mpp ) CALL mpp_sum( zq2 ) 991 zinvt = zq1 + zq2 992 zdelta = zinvt - cycletot2(6) 993 zratio = 1.0e2 * zdelta / cycletot2(6) 994 ! 995 IF ( lwp ) WRITE(numout,9010) 'oxygen', zq1, zq2, zinvt, & 996 cycletot2(6), zdelta, zratio 997 IF ( lwp ) WRITE(numout,*) 998 ! 999 !--------------------------------------------------------------------- 1000 zq1 = loc_vol 1001 zq2 = loc_area 1002 IF( lk_mpp ) CALL mpp_sum( zq1 ) 1003 IF( lk_mpp ) CALL mpp_sum( zq2 ) 1004 IF ( lwp ) WRITE(numout,*) 1005 IF ( lwp ) WRITE(numout,*) ' check : cvol : ', zq1 1006 IF ( lwp ) WRITE(numout,*) ' check : carea : ', zq2 1007 IF ( lwp ) WRITE(numout,*) 1008 IF ( lwp ) CALL flush(numout) 1009 1010 !! 1011 !! 1012 9010 FORMAT(' element:',a10, & 1013 ' 3d sum:',e18.10,' 2d sum:',e18.10, & 1014 ' total:',e18.10,' initial:',e18.10, & 1015 ' delta:',e18.10,' %:',e18.10) 1016 ! 1017 END SUBROUTINE trc_rst_conserve 1018 #endif 1019 699 1020 SUBROUTINE trc_rst_tra_stat 700 1021 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
