Changeset 15600 for NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate
- Timestamp:
- 2021-12-15T12:49:49+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90
r15590 r15600 752 752 753 753 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: amp_u2d,phi_u2d, amp_v2d,phi_v2d ! arrays for output 754 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:,:) :: amp_u3d,phi_u3d, amp_v3d,phi_v3d ! arrays for output 754 755 755 756 REAL(wp) :: tmp_u_amp ,tmp_v_amp ,tmp_u_phi ,tmp_v_phi … … 758 759 REAL(wp) :: tmpreal 759 760 760 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: tmp_u_amp_mat,tmp_v_amp_mat,tmp_u_phi_mat,tmp_v_phi_mat 761 ! REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: a_u_mat,b_u_mat,a_v_mat,b_v_mat,qmax_mat,qmin_mat,ecc_mat 762 ! REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: thetamax_mat,thetamin_mat,Qc_mat,Qac_mat,gc_mat,gac_mat 763 ! REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: Phi_Ua_mat,dir_Ua_mat,polarity_mat 764 765 766 767 ! IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN 768 ! ALLOCATE( amp_u2d(jh,jpi,jpj),amp_v2d(jh,jpi,jpj),phi_u2d(jh,jpi,jpj),phi_v2d(jh,jpi,jpj) ) 769 770 771 ! ALLOCATE(tmp_u_amp_mat(jpi,jpj),tmp_v_amp_mat(jpi,jpj),tmp_u_phi_mat(jpi,jpj),tmp_v_phi_mat(jpi,jpj)) 772 !! ALLOCATE(a_u_mat(jpi,jpj),b_u_mat(jpi,jpj),a_v_mat(jpi,jpj),b_v_mat(jpi,jpj)) 773 !! ALLOCATE(qmax_mat(jpi,jpj),qmin_mat(jpi,jpj),ecc_mat(jpi,jpj)) 774 !! ALLOCATE(thetamax_mat(jpi,jpj),thetamin_mat(jpi,jpj),Qc_mat(jpi,jpj),Qac_mat(jpi,jpj)) 775 !! ALLOCATE(gc_mat(jpi,jpj),gac_mat(jpi,jpj),Phi_Ua_mat(jpi,jpj),dir_Ua_mat(jpi,jpj),polarity_mat(jpi,jpj)) 776 777 ! endif 761 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: tmp_u_amp_2d_mat,tmp_v_amp_2d_mat,tmp_u_phi_2d_mat,tmp_v_phi_2d_mat 762 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: a_u_2d_mat,b_u_2d_mat,a_v_2d_mat,b_v_2d_mat 763 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: qmax_2d_mat,qmin_2d_mat,ecc_2d_mat 764 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: thetamax_2d_mat,thetamin_2d_mat,Qc_2d_mat,Qac_2d_mat 765 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: gc_2d_mat,gac_2d_mat,Phi_Ua_2d_mat,dir_Ua_2d_mat 766 REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: polarity_2d_mat 767 768 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: tmp_u_amp_3d_mat,tmp_v_amp_3d_mat,tmp_u_phi_3d_mat,tmp_v_phi_3d_mat 769 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: a_u_3d_mat,b_u_3d_mat,a_v_3d_mat,b_v_3d_mat 770 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: qmax_3d_mat,qmin_3d_mat,ecc_3d_mat 771 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: thetamax_3d_mat,thetamin_3d_mat,Qc_3d_mat,Qac_3d_mat 772 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: gc_3d_mat,gac_3d_mat,Phi_Ua_3d_mat,dir_Ua_3d_mat 773 REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: polarity_3d_mat 774 775 776 777 778 IF (ln_diaharm_postproc_vel) THEN 779 IF (ln_ana_uvbar) THEN 780 ALLOCATE( amp_u2d(nb_ana,jpi,jpj), amp_v2d(nb_ana,jpi,jpj), phi_u2d(nb_ana,jpi,jpj), phi_v2d(nb_ana,jpi,jpj) ) 781 782 783 ALLOCATE(tmp_u_amp_2d_mat(jpi,jpj),tmp_v_amp_2d_mat(jpi,jpj),tmp_u_phi_2d_mat(jpi,jpj),tmp_v_phi_2d_mat(jpi,jpj)) 784 ALLOCATE(a_u_2d_mat(jpi,jpj),b_u_2d_mat(jpi,jpj),a_v_2d_mat(jpi,jpj),b_v_2d_mat(jpi,jpj)) 785 ALLOCATE(qmax_2d_mat(jpi,jpj),qmin_2d_mat(jpi,jpj),ecc_2d_mat(jpi,jpj)) 786 ALLOCATE(thetamax_2d_mat(jpi,jpj),thetamin_2d_mat(jpi,jpj),Qc_2d_mat(jpi,jpj),Qac_2d_mat(jpi,jpj)) 787 ALLOCATE(gc_2d_mat(jpi,jpj),gac_2d_mat(jpi,jpj),Phi_Ua_2d_mat(jpi,jpj),dir_Ua_2d_mat(jpi,jpj)) 788 ALLOCATE(polarity_2d_mat(jpi,jpj)) 789 790 ENDIF 791 792 793 IF (ln_ana_uv3d) THEN 794 ALLOCATE( amp_u3d(nb_ana,jpi,jpj,jpk), amp_v3d(nb_ana,jpi,jpj,jpk), phi_u3d(nb_ana,jpi,jpj,jpk), phi_v3d(nb_ana,jpi,jpj,jpk) ) 795 796 797 ALLOCATE(tmp_u_amp_3d_mat(jpi,jpj,jpk),tmp_v_amp_3d_mat(jpi,jpj,jpk),tmp_u_phi_3d_mat(jpi,jpj,jpk),tmp_v_phi_3d_mat(jpi,jpj,jpk)) 798 ALLOCATE(a_u_3d_mat(jpi,jpj,jpk),b_u_3d_mat(jpi,jpj,jpk),a_v_3d_mat(jpi,jpj,jpk),b_v_3d_mat(jpi,jpj,jpk)) 799 ALLOCATE(qmax_3d_mat(jpi,jpj,jpk),qmin_3d_mat(jpi,jpj,jpk),ecc_3d_mat(jpi,jpj,jpk)) 800 ALLOCATE(thetamax_3d_mat(jpi,jpj,jpk),thetamin_3d_mat(jpi,jpj,jpk),Qc_3d_mat(jpi,jpj,jpk),Qac_3d_mat(jpi,jpj,jpk)) 801 ALLOCATE(gc_3d_mat(jpi,jpj,jpk),gac_3d_mat(jpi,jpj,jpk),Phi_Ua_3d_mat(jpi,jpj,jpk),dir_Ua_3d_mat(jpi,jpj,jpk)) 802 ALLOCATE(polarity_3d_mat(jpi,jpj,jpk)) 803 804 ENDIF 805 ENDIF 778 806 779 807 do jgrid=1,nvar_2d … … 810 838 ! NETCDF OUTPUT 811 839 suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) ) 812 IF(lwp) WRITE(numout,*) " harm_ana_out", suffix840 IF(lwp) WRITE(numout,*) "diaharm_fast", suffix 813 841 814 842 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 815 843 IF( iom_use(TRIM(tmp_name)) ) THEN 816 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)817 IF(lwp) WRITE(numout,*) " harm_ana_out names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"844 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D) 845 IF(lwp) WRITE(numout,*) "diaharm_fast names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr" 818 846 CALL iom_put( TRIM(tmp_name), h_out2D(:,:) ) 819 847 ELSE 820 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)848 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 821 849 ENDIF 822 850 823 851 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 824 852 IF( iom_use(TRIM(tmp_name)) ) THEN 825 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)853 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D) 826 854 CALL iom_put( TRIM(tmp_name), g_out2D(:,:) ) 827 855 ELSE 828 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)856 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 829 857 ENDIF 830 858 831 859 832 860 833 ! IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN 834 835 ! !IF (m_posi_2d(jgrid) == 2) THEN 836 ! IF (TRIM(suffix) == TRIM('u2d')) THEN 837 ! if (lwp) WRITE(numout,*) "harm_ana_out ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u2d '//TRIM(suffix) 838 ! amp_u2d(jh,:,:) = h_out2D(:,:) 839 ! phi_u2d(jh,:,:) = rpi*g_out2D(:,:)/180.0 840 ! ENDIF 841 842 ! !IF (m_posi_2d(jgrid) == 3) THEN 843 ! IF (TRIM(suffix) == TRIM('v2d')) THEN 844 ! if (lwp) WRITE(numout,*) "harm_ana_out ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v2d '//TRIM(suffix) 845 ! amp_v2d(jh,:,:) = h_out2D(:,:) 846 ! phi_v2d(jh,:,:) = rpi*g_out2D(:,:)/180.0 847 ! ENDIF 848 ! ENDIF 849 850 ! CALL FLUSH(numout) 861 IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN 862 863 !IF (m_posi_2d(jgrid) == 2) THEN 864 IF (TRIM(suffix) == TRIM('u2d')) THEN 865 if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u2d '//TRIM(suffix) 866 do jj=1,nlcj 867 do ji=1,nlci 868 if (ssumask(ji,jj) == 1) THEN 869 amp_u2d(jh,ji,jj) = h_out2D(ji,jj) 870 phi_u2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0 871 else 872 amp_u2d(jh,ji,jj) = 0. 873 phi_u2d(jh,ji,jj) = 0. 874 ENDIF 875 enddo 876 enddo 877 ENDIF 878 879 !IF (m_posi_2d(jgrid) == 3) THEN 880 IF (TRIM(suffix) == TRIM('v2d')) THEN 881 if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v2d '//TRIM(suffix) 882 do jj=1,nlcj 883 do ji=1,nlci 884 if (ssvmask(ji,jj) == 1) THEN 885 amp_v2d(jh,ji,jj) = h_out2D(ji,jj) 886 phi_v2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0 887 else 888 amp_v2d(jh,ji,jj) = 0 889 phi_v2d(jh,ji,jj) = 0 890 ENDIF 891 enddo 892 enddo 893 ENDIF 894 ENDIF 895 896 CALL FLUSH(numout) 851 897 852 898 … … 856 902 tmp_name='TA_'//TRIM(suffix)//'_off' 857 903 IF( iom_use(TRIM(tmp_name)) ) THEN 858 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name)904 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 859 905 CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid)) 860 906 ELSE 861 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)907 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 862 908 ENDIF 863 909 … … 903 949 ! NETCDF OUTPUT 904 950 suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) ) 905 IF(lwp) WRITE(numout,*) " harm_ana_out", suffix951 IF(lwp) WRITE(numout,*) "diaharm_fast", suffix 906 952 907 953 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 908 954 IF( iom_use(TRIM(tmp_name)) ) THEN 909 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)955 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D) 910 956 CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) ) 911 957 ELSE 912 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)958 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 913 959 ENDIF 914 960 915 961 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 916 962 IF( iom_use(TRIM(tmp_name)) ) THEN 917 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)963 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D) 918 964 CALL iom_put(tmp_name, g_out3D(:,:,:) ) 919 965 ELSE 920 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)966 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 921 967 ENDIF 922 968 … … 926 972 tmp_name='TA_'//TRIM(suffix)//'_off' 927 973 IF( iom_use(TRIM(tmp_name)) ) THEN 928 IF(lwp) WRITE(numout,*) " harm_ana_out: iom_put: ",TRIM(tmp_name)974 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 929 975 CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid)) 930 976 ELSE 931 IF(lwp) WRITE(numout,*) " harm_ana_out: not requested: ",TRIM(tmp_name)977 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 932 978 ENDIF 933 979 980 981 982 983 984 985 IF (ln_diaharm_postproc_vel .AND. ln_ana_uv3d) THEN 986 987 !IF (m_posi_2d(jgrid) == 2) THEN 988 IF (TRIM(suffix) == TRIM('u3d')) THEN 989 if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u3d '//TRIM(suffix) 990 DO jk=1,jpkm1 991 do jj=1,nlcj 992 do ji=1,nlci 993 if (umask(ji,jj,jk) == 1) THEN 994 amp_u3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) 995 phi_u3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0 996 else 997 amp_u3d(jh,ji,jj,jk) = 0. 998 phi_u3d(jh,ji,jj,jk) = 0. 999 ENDIF 1000 enddo 1001 enddo 1002 enddo 1003 ENDIF 1004 1005 !IF (m_posi_2d(jgrid) == 3) THEN 1006 IF (TRIM(suffix) == TRIM('v3d')) THEN 1007 if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v3d '//TRIM(suffix) 1008 DO jk=1,jpkm1 1009 do jj=1,nlcj 1010 do ji=1,nlci 1011 if (vmask(ji,jj,jk) == 1) THEN 1012 amp_v3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) 1013 phi_v3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0 1014 else 1015 amp_v3d(jh,ji,jj,jk) = 0 1016 phi_v3d(jh,ji,jj,jk) = 0 1017 ENDIF 1018 enddo 1019 enddo 1020 enddo 1021 ENDIF 1022 ENDIF 1023 1024 CALL FLUSH(numout) 1025 934 1026 enddo ! jgrid 935 1027 936 1028 CALL FLUSH(numout) 937 1029 938 ! IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN 939 ! IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess barotropic velocity tidal parameters" 940 ! CALL FLUSH(numout) 941 ! DO jh=1,nb_ana 942 943 944 ! tmp_u_amp_mat(:,:) = 0. 945 ! tmp_v_amp_mat(:,:) = 0. 946 ! tmp_u_phi_mat(:,:) = 0. 947 ! tmp_v_phi_mat(:,:) = 0. 948 949 !! a_u_mat(:,:) = 0. 950 !! b_u_mat(:,:) = 0. 951 !! a_v_mat(:,:) = 0. 952 !! b_v_mat(:,:) = 0. 953 954 !! qmax_mat(:,:) = 0. 955 !! qmin_mat(:,:) = 0. 956 957 !! ecc_mat(:,:) = 0 958 !! thetamax_mat(:,:) =0. 959 !! thetamin_mat(:,:) = 0. 960 961 !! Qc_mat(:,:) = 0. 962 !! Qac_mat(:,:) = 0. 963 !! gc_mat(:,:) = 0. 964 !! gac_mat(:,:) = 0. 965 966 !! Phi_Ua_mat(:,:) = 0. 967 !! dir_Ua_mat(:,:) = 0. 968 !! polarity_mat(:,:) = 0. 969 970 971 !! DO jj = 2, nlcj - 1 972 !! DO ji = 2, nlci - 1 973 974 !! do jj=2,nlcj 975 !! do ji=2,nlci 976 ! !IF ((ssumask(ji,jj) + ssumask(ji-1,jj)) == 0 ) CYCLE 977 ! !IF ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) == 0 ) CYCLE 978 979 !! IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN 980 !! tmp_u_amp = ((amp_u2d(jh,ji,jj)*ssumask(ji,jj)) + (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) 981 !! tmp_v_amp = ((amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) 982 !! ! WORK ON THE WRAP AROUND 983 !! tmp_u_phi = ((phi_u2d(jh,ji,jj)*ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) 984 !! tmp_v_phi = ((phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) 985 986 ! do jj=1,nlcj 987 ! do ji=1,nlci 988 989 !! tmp_u_amp = ((amp_u2d(jh,ji,jj)) + (amp_u2d(jh,ji-1,jj)))/(2.) 990 !! tmp_v_amp = ((amp_v2d(jh,ji,jj)) + (amp_v2d(jh,ji,jj-1)))/(2.) 991 !! ! WORK ON THE WRAP AROUND 992 !! tmp_u_phi = ((phi_u2d(jh,ji,jj)) + (phi_u2d(jh,ji-1,jj)))/(2.) 993 !! tmp_v_phi = ((phi_v2d(jh,ji,jj)) + (phi_v2d(jh,ji,jj-1)))/(2.) 994 995 996 997 ! tmp_u_amp = (amp_u2d(jh,ji,jj)) 998 ! tmp_v_amp = (amp_v2d(jh,ji,jj)) 999 ! ! WORK ON THE WRAP AROUND 1000 ! tmp_u_phi = (phi_u2d(jh,ji,jj)) 1001 ! tmp_v_phi = (phi_v2d(jh,ji,jj)) 1002 1003 1004 1005 !! a_u = tmp_U_amp * cos(tmp_U_phi) 1006 !! b_u = tmp_U_amp * sin(tmp_U_phi) 1007 !! a_v = tmp_V_amp * cos(tmp_V_phi) 1008 !! b_v = tmp_V_amp * sin(tmp_V_phi) 1009 1010 !! twodelta = atan2( (tmp_V_amp**2 * sin( 2*(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp**2 + tmp_V_amp**2 * cos( 2*(tmp_U_phi - tmp_V_phi) ) ) ) 1011 !! delta = twodelta/2. 1012 1013 !! !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1014 1015 !! tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 1016 !! if (tmpreal < 0) CYCLE 1017 !! alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1018 !! if (alpha2 < 0) CYCLE 1019 !! alpha= sqrt( alpha2 ) 1020 1021 1022 !! !major and minor axis of the ellipse 1023 !! !qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) 1024 !! !tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1025 !! !qmin = 0 1026 !! !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1027 1028 !! tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1029 !! if (tmpreal < 0) CYCLE 1030 !! qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1031 1032 !! !eccentricity of ellipse 1033 !! ecc = (qmax - qmin)/(qmax + qmin) 1034 !! ! Angle of major and minor ellipse 1035 !! thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) 1036 !! thetamin = thetamax + rpi/2. 1037 1038 1039 1040 !! ! Rotary current components: Pugh A3.10 1041 !! ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors 1042 !! ! so Qc = clockwise = anticyclonic = negative 1043 !! ! and Qac = anticlockwise = cyclonic = negative 1044 1045 !! tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1046 !! if (tmpreal < 0) CYCLE 1047 !! Qc = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1048 1049 !! tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1050 !! if (tmpreal < 0) CYCLE 1051 !! Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1052 1053 1054 !! gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1055 !! gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1056 1057 !! !Pugh A3.2 1058 !! Phi_Ua = -0.5*(gac - gc) 1059 !! dir_Ua = 0.5*(gac + gc) ! positive from x axis 1060 !! polarity = (Qac - Qc)/qmax 1061 1062 1063 1064 ! tmp_u_amp_mat(ji,jj) = tmp_u_amp 1065 ! tmp_v_amp_mat(ji,jj) = tmp_v_amp 1066 ! tmp_u_phi_mat(ji,jj) = tmp_u_phi 1067 ! tmp_v_phi_mat(ji,jj) = tmp_v_phi 1068 1069 1070 !! a_u_mat(ji,jj) = a_u 1071 !! b_u_mat(ji,jj) = b_u 1072 !! a_v_mat(ji,jj) = a_v 1073 !! b_v_mat(ji,jj) = b_v 1074 1075 !! qmax_mat(ji,jj) = qmax 1076 !! qmin_mat(ji,jj) = qmin 1077 1078 !! ecc_mat(ji,jj) = ecc 1079 !! thetamax_mat(ji,jj) = thetamax 1080 !! thetamin_mat(ji,jj) = thetamin 1081 1082 !! Qc_mat(ji,jj) = Qc 1083 !! Qac_mat(ji,jj) = Qac 1084 !! gc_mat(ji,jj) = gc 1085 !! gac_mat(ji,jj) = gac 1086 1087 !! Phi_Ua_mat(ji,jj) = Phi_Ua 1088 !! dir_Ua_mat(ji,jj) = dir_Ua 1089 !! polarity_mat(ji,jj) = polarity 1090 1091 !! ENDIF 1092 ! END DO 1093 ! END DO 1094 1095 1096 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar' 1097 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1098 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1099 !! CALL iom_put( TRIM(tmp_name), tmp_u_amp_mat(:,:)) 1100 !! ENDIF 1101 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar' 1102 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1103 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1104 !! CALL iom_put( TRIM(tmp_name), tmp_v_amp_mat(:,:)) 1105 !! ENDIF 1106 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar' 1107 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1108 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1109 !! CALL iom_put( TRIM(tmp_name), tmp_u_phi_mat(:,:)) 1110 !! ENDIF 1111 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar' 1112 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1113 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1114 !! CALL iom_put( TRIM(tmp_name), tmp_v_phi_mat(:,:)) 1115 !! ENDIF 1116 1117 1118 1119 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar' 1120 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1121 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1122 !! CALL iom_put( TRIM(tmp_name), a_u_mat(:,:)) 1123 !! ENDIF 1124 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar' 1125 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1126 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1127 !! CALL iom_put( TRIM(tmp_name), a_v_mat(:,:)) 1128 !! ENDIF 1129 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar' 1130 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1131 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1132 !! CALL iom_put( TRIM(tmp_name), b_u_mat(:,:)) 1133 !! ENDIF 1134 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar' 1135 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1136 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1137 !! CALL iom_put( TRIM(tmp_name), b_v_mat(:,:)) 1138 !! ENDIF 1139 1140 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar' 1141 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1142 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1143 !! CALL iom_put( TRIM(tmp_name), qmax_mat(:,:)) 1144 !! ENDIF 1145 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar' 1146 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1147 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1148 !! CALL iom_put( TRIM(tmp_name), qmin_mat(:,:)) 1149 !! ENDIF 1150 1151 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar' 1152 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1153 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1154 !! CALL iom_put( TRIM(tmp_name), ecc_mat(:,:)) 1155 !! ENDIF 1156 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar' 1157 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1158 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1159 !! CALL iom_put( TRIM(tmp_name), thetamax_mat(:,:)) 1160 !! ENDIF 1161 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar' 1162 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1163 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1164 !! CALL iom_put( TRIM(tmp_name), thetamin_mat(:,:)) 1165 !! ENDIF 1166 1167 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar' 1168 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1169 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1170 !! CALL iom_put( TRIM(tmp_name), Qc_mat(:,:)) 1171 !! ENDIF 1172 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar' 1173 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1174 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1175 !! CALL iom_put( TRIM(tmp_name), Qac_mat(:,:)) 1176 !! ENDIF 1177 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar' 1178 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1179 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1180 !! CALL iom_put( TRIM(tmp_name), gc_mat(:,:)) 1181 !! ENDIF 1182 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar' 1183 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1184 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1185 !! CALL iom_put( TRIM(tmp_name), gac_mat(:,:)) 1186 !! ENDIF 1187 1188 1189 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar' 1190 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1191 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1192 !! CALL iom_put( TRIM(tmp_name), Phi_Ua_mat(:,:)) 1193 !! ENDIF 1194 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar' 1195 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1196 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1197 !! CALL iom_put( TRIM(tmp_name), dir_Ua_mat(:,:)) 1198 !! ENDIF 1199 !! tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar' 1200 !! IF( iom_use(TRIM(tmp_name)) ) THEN 1201 !! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name) 1202 !! CALL iom_put( TRIM(tmp_name), polarity_mat(:,:)) 1203 !! ENDIF 1204 1205 ! tmp_u_amp_mat(:,:) = 0. 1206 ! tmp_v_amp_mat(:,:) = 0. 1207 ! tmp_u_phi_mat(:,:) = 0. 1208 ! tmp_v_phi_mat(:,:) = 0. 1209 1210 !! a_u_mat(:,:) = 0. 1211 !! b_u_mat(:,:) = 0. 1212 !! a_v_mat(:,:) = 0. 1213 !! b_v_mat(:,:) = 0. 1214 1215 !! qmax_mat(:,:) = 0. 1216 !! qmin_mat(:,:) = 0. 1217 1218 !! ecc_mat(:,:) = 0 1219 !! thetamax_mat(:,:) =0. 1220 !! thetamin_mat(:,:) = 0. 1221 1222 !! Qc_mat(:,:) = 0. 1223 !! Qac_mat(:,:) = 0. 1224 !! gc_mat(:,:) = 0. 1225 !! gac_mat(:,:) = 0. 1226 1227 !! Phi_Ua_mat(:,:) = 0. 1228 !! dir_Ua_mat(:,:) = 0. 1229 !! polarity_mat(:,:) = 0. 1230 1231 1232 ! END DO 1233 ! IF(lwp) WRITE(numout,*) "diaharm_fast: Finshed postprocessing 2d velocity tidal parameters" 1234 ! ENDIF 1235 1236 ! CALL FLUSH(numout) 1237 1238 ! IF (ln_diaharm_postproc_vel .AND. ln_ana_uv3d) THEN 1239 ! IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess 3d velocity tidal parameters" 1240 ! ENDIF 1241 1242 1243 ! CALL FLUSH(numout) 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 IF (ln_diaharm_postproc_vel ) THEN 1047 IF ( ln_ana_uvbar) THEN 1048 IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess barotropic velocity tidal parameters" 1049 CALL FLUSH(numout) 1050 DO jh=1,nb_ana 1051 1052 1053 tmp_u_amp_2d_mat(:,:) = 0. 1054 tmp_v_amp_2d_mat(:,:) = 0. 1055 tmp_u_phi_2d_mat(:,:) = 0. 1056 tmp_v_phi_2d_mat(:,:) = 0. 1057 1058 a_u_2d_mat(:,:) = 0. 1059 b_u_2d_mat(:,:) = 0. 1060 a_v_2d_mat(:,:) = 0. 1061 b_v_2d_mat(:,:) = 0. 1062 1063 qmax_2d_mat(:,:) = 0. 1064 qmin_2d_mat(:,:) = 0. 1065 1066 ecc_2d_mat(:,:) = 0. 1067 thetamax_2d_mat(:,:) =0. 1068 thetamin_2d_mat(:,:) = 0. 1069 1070 Qc_2d_mat(:,:) = 0. 1071 Qac_2d_mat(:,:) = 0. 1072 gc_2d_mat(:,:) = 0. 1073 gac_2d_mat(:,:) = 0. 1074 1075 Phi_Ua_2d_mat(:,:) = 0. 1076 dir_Ua_2d_mat(:,:) = 0. 1077 polarity_2d_mat(:,:) = 0. 1078 1079 1080 DO jj = 2, nlcj !- 1 1081 DO ji = 2, nlci ! - 1 1082 1083 ! do jj=2,nlcj 1084 ! do ji=2,nlci 1085 !IF ((ssumask(ji,jj) + ssumask(ji-1,jj)) == 0 ) CYCLE 1086 !IF ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) == 0 ) CYCLE 1087 1088 IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN 1089 tmp_u_amp = ((amp_u2d(jh,ji,jj)*ssumask(ji,jj)) + (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) 1090 tmp_v_amp = ((amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) 1091 ! WORK ON THE WRAP AROUND 1092 !tmp_u_phi = ((phi_u2d(jh,ji,jj)*ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) 1093 !tmp_v_phi = ((phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) 1094 1095 if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then 1096 !tmp_u_phi = ((phi_u2d(jh,ji,jj)*ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) 1097 tmp_u_phi = atan2((sin(phi_u2d(jh,ji,jj)) + sin(phi_u2d(jh,ji-1,jj))),(cos(phi_u2d(jh,ji,jj)) + cos(phi_u2d(jh,ji-1,jj)))) 1098 else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then 1099 tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj)) 1100 else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then 1101 tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) 1102 else 1103 cycle 1104 end if 1105 1106 1107 if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then 1108 !tmp_v_phi = ((phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) 1109 tmp_v_phi = atan2((sin(phi_v2d(jh,ji,jj)) + sin(phi_v2d(jh,ji,jj-1))),(cos(phi_v2d(jh,ji,jj)) + cos(phi_v2d(jh,ji,jj-1)))) 1110 else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then 1111 tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) 1112 else if ( (ssvmask(ji,jj) == 0) .AND. (ssvmask(ji,jj-1) == 1)) then 1113 tmp_v_phi = (phi_v2d(ji,jj-1)*ssvmask(ji,jj-1)) 1114 else 1115 cycle 1116 end if 1117 1118 1119 ! do jj=1,nlcj 1120 ! do ji=1,nlci 1121 1122 ! tmp_u_amp = ((amp_u2d(jh,ji,jj)) + (amp_u2d(jh,ji-1,jj)))/(2.) 1123 ! tmp_v_amp = ((amp_v2d(jh,ji,jj)) + (amp_v2d(jh,ji,jj-1)))/(2.) 1124 ! ! WORK ON THE WRAP AROUND 1125 ! tmp_u_phi = ((phi_u2d(jh,ji,jj)) + (phi_u2d(jh,ji-1,jj)))/(2.) 1126 ! tmp_v_phi = ((phi_v2d(jh,ji,jj)) + (phi_v2d(jh,ji,jj-1)))/(2.) 1127 1128 1129 1130 ! tmp_u_amp = (amp_u2d(jh,ji,jj)) 1131 ! tmp_v_amp = (amp_v2d(jh,ji,jj)) 1132 ! ! WORK ON THE WRAP AROUND 1133 ! tmp_u_phi = (phi_u2d(jh,ji,jj)) 1134 ! tmp_v_phi = (phi_v2d(jh,ji,jj)) 1135 1136 1137 1138 ! a_u = tmp_U_amp * cos(tmp_U_phi) 1139 ! b_u = tmp_U_amp * sin(tmp_U_phi) 1140 ! a_v = tmp_V_amp * cos(tmp_V_phi) 1141 ! b_v = tmp_V_amp * sin(tmp_V_phi) 1142 1143 ! twodelta = atan2( (tmp_V_amp**2 * sin( 2*(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp**2 + tmp_V_amp**2 * cos( 2*(tmp_U_phi - tmp_V_phi) ) ) ) 1144 ! delta = twodelta/2. 1145 1146 ! !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1147 1148 ! tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 1149 ! if (tmpreal < 0) CYCLE 1150 ! alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1151 ! if (alpha2 < 0) CYCLE 1152 ! alpha= sqrt( alpha2 ) 1153 1154 1155 ! !major and minor axis of the ellipse 1156 ! qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) 1157 ! !tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1158 ! !qmin = 0 1159 ! !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1160 1161 ! tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1162 ! if (tmpreal < 0) CYCLE 1163 ! qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1164 1165 ! !eccentricity of ellipse 1166 ! tmpreal = (qmax + qmin) 1167 ! if (tmpreal < 0) CYCLE 1168 ! ecc = (qmax - qmin)/(qmax + qmin) 1169 ! ! Angle of major and minor ellipse 1170 ! thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) 1171 ! thetamin = thetamax + rpi/2. 1172 1173 1174 1175 ! ! Rotary current components: Pugh A3.10 1176 ! ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors 1177 ! ! so Qc = clockwise = anticyclonic = negative 1178 ! ! and Qac = anticlockwise = cyclonic = negative 1179 1180 ! tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1181 ! if (tmpreal < 0) CYCLE 1182 ! Qc = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1183 1184 ! tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1185 ! if (tmpreal < 0) CYCLE 1186 ! Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1187 1188 1189 ! gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1190 ! gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1191 1192 ! !Pugh A3.2 1193 ! Phi_Ua = -0.5*(gac - gc) 1194 ! dir_Ua = 0.5*(gac + gc) ! positive from x axis 1195 1196 ! tmpreal = qmax 1197 ! if (tmpreal < 0) CYCLE 1198 ! polarity = (Qac - Qc)/qmax 1199 1200 a_u = tmp_U_amp * cos(tmp_U_phi) 1201 b_u = tmp_U_amp * sin(tmp_U_phi) 1202 a_v = tmp_V_amp * cos(tmp_V_phi) 1203 b_v = tmp_V_amp * sin(tmp_V_phi) 1204 1205 twodelta = atan2( (tmp_V_amp**2 * sin( 2*(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp**2 + tmp_V_amp**2 * cos( 2*(tmp_U_phi - tmp_V_phi) ) ) ) 1206 delta = twodelta/2. 1207 1208 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1209 1210 tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 1211 if (tmpreal < 0) tmpreal = 0 !CYCLE 1212 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1213 alpha2 = sqrt( tmpreal ) 1214 if (alpha2 < 0) alpha2 = 0 !CYCLE 1215 alpha= sqrt( alpha2 ) 1216 1217 1218 !major and minor axis of the ellipse 1219 qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) 1220 !tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1221 !qmin = 0 1222 !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1223 1224 tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1225 if (tmpreal < 0) tmpreal = 0 !CYCLE 1226 !qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1227 qmin = sqrt( tmpreal ) ! but always positive. 1228 1229 !eccentricity of ellipse 1230 1231 tmpreal = (qmax + qmin) 1232 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE 1233 !ecc = (qmax - qmin)/(qmax + qmin) 1234 ecc = (qmax - qmin)/(tmpreal) 1235 ! Angle of major and minor ellipse 1236 thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) 1237 thetamin = thetamax + rpi/2. 1238 1239 1240 1241 ! Rotary current components: Pugh A3.10 1242 ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors 1243 ! so Qc = clockwise = anticyclonic = negative 1244 ! and Qac = anticlockwise = cyclonic = negative 1245 1246 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1247 if (tmpreal < 0) tmpreal = 0! CYCLE 1248 !Qc = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1249 Qc = 0.5*sqrt( tmpreal ) 1250 1251 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1252 if (tmpreal < 0) tmpreal = 0! CYCLE 1253 !Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1254 Qac = 0.5*sqrt( tmpreal ) 1255 1256 1257 gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1258 gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1259 1260 !Pugh A3.2 1261 Phi_Ua = -0.5*(gac - gc) 1262 dir_Ua = 0.5*(gac + gc) ! positive from x axis 1263 1264 tmpreal = qmax 1265 !if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 1266 if (tmpreal == 0) then 1267 polarity = 0 1268 else 1269 polarity = (Qac - Qc)/qmax 1270 endif 1271 1272 1273 tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp 1274 tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp 1275 tmp_u_phi_2d_mat(ji,jj) = tmp_u_phi 1276 tmp_v_phi_2d_mat(ji,jj) = tmp_v_phi 1277 1278 1279 a_u_2d_mat(ji,jj) = a_u 1280 b_u_2d_mat(ji,jj) = b_u 1281 a_v_2d_mat(ji,jj) = a_v 1282 b_v_2d_mat(ji,jj) = b_v 1283 1284 qmax_2d_mat(ji,jj) = qmax 1285 qmin_2d_mat(ji,jj) = qmin 1286 1287 ecc_2d_mat(ji,jj) = ecc 1288 thetamax_2d_mat(ji,jj) = thetamax 1289 thetamin_2d_mat(ji,jj) = thetamin 1290 1291 Qc_2d_mat(ji,jj) = Qc 1292 Qac_2d_mat(ji,jj) = Qac 1293 gc_2d_mat(ji,jj) = gc 1294 gac_2d_mat(ji,jj) = gac 1295 1296 Phi_Ua_2d_mat(ji,jj) = Phi_Ua 1297 dir_Ua_2d_mat(ji,jj) = dir_Ua 1298 polarity_2d_mat(ji,jj) = polarity 1299 1300 ENDIF 1301 END DO 1302 END DO 1303 1304 1305 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar' 1306 IF( iom_use(TRIM(tmp_name)) ) THEN 1307 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1308 CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:)) 1309 ENDIF 1310 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar' 1311 IF( iom_use(TRIM(tmp_name)) ) THEN 1312 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1313 CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:)) 1314 ENDIF 1315 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar' 1316 IF( iom_use(TRIM(tmp_name)) ) THEN 1317 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1318 CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:)) 1319 ENDIF 1320 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar' 1321 IF( iom_use(TRIM(tmp_name)) ) THEN 1322 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1323 CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:)) 1324 ENDIF 1325 1326 1327 1328 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar' 1329 IF( iom_use(TRIM(tmp_name)) ) THEN 1330 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1331 CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:)) 1332 ENDIF 1333 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar' 1334 IF( iom_use(TRIM(tmp_name)) ) THEN 1335 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1336 CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:)) 1337 ENDIF 1338 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar' 1339 IF( iom_use(TRIM(tmp_name)) ) THEN 1340 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1341 CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:)) 1342 ENDIF 1343 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar' 1344 IF( iom_use(TRIM(tmp_name)) ) THEN 1345 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1346 CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:)) 1347 ENDIF 1348 1349 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar' 1350 IF( iom_use(TRIM(tmp_name)) ) THEN 1351 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1352 CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:)) 1353 ENDIF 1354 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar' 1355 IF( iom_use(TRIM(tmp_name)) ) THEN 1356 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1357 CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:)) 1358 ENDIF 1359 1360 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar' 1361 IF( iom_use(TRIM(tmp_name)) ) THEN 1362 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1363 CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:)) 1364 ENDIF 1365 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar' 1366 IF( iom_use(TRIM(tmp_name)) ) THEN 1367 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1368 CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:)) 1369 ENDIF 1370 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar' 1371 IF( iom_use(TRIM(tmp_name)) ) THEN 1372 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1373 CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:)) 1374 ENDIF 1375 1376 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar' 1377 IF( iom_use(TRIM(tmp_name)) ) THEN 1378 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1379 CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:)) 1380 ENDIF 1381 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar' 1382 IF( iom_use(TRIM(tmp_name)) ) THEN 1383 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1384 CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:)) 1385 ENDIF 1386 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar' 1387 IF( iom_use(TRIM(tmp_name)) ) THEN 1388 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1389 CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:)) 1390 ENDIF 1391 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar' 1392 IF( iom_use(TRIM(tmp_name)) ) THEN 1393 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1394 CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:)) 1395 ENDIF 1396 1397 1398 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar' 1399 IF( iom_use(TRIM(tmp_name)) ) THEN 1400 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1401 CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:)) 1402 ENDIF 1403 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar' 1404 IF( iom_use(TRIM(tmp_name)) ) THEN 1405 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1406 CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:)) 1407 ENDIF 1408 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar' 1409 IF( iom_use(TRIM(tmp_name)) ) THEN 1410 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1411 CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:)) 1412 ENDIF 1413 1414 tmp_u_amp_2d_mat(:,:) = 0. 1415 tmp_v_amp_2d_mat(:,:) = 0. 1416 tmp_u_phi_2d_mat(:,:) = 0. 1417 tmp_v_phi_2d_mat(:,:) = 0. 1418 1419 a_u_2d_mat(:,:) = 0. 1420 b_u_2d_mat(:,:) = 0. 1421 a_v_2d_mat(:,:) = 0. 1422 b_v_2d_mat(:,:) = 0. 1423 1424 qmax_2d_mat(:,:) = 0. 1425 qmin_2d_mat(:,:) = 0. 1426 1427 ecc_2d_mat(:,:) = 0. 1428 thetamax_2d_mat(:,:) =0. 1429 thetamin_2d_mat(:,:) = 0. 1430 1431 Qc_2d_mat(:,:) = 0. 1432 Qac_2d_mat(:,:) = 0. 1433 gc_2d_mat(:,:) = 0. 1434 gac_2d_mat(:,:) = 0. 1435 1436 Phi_Ua_2d_mat(:,:) = 0. 1437 dir_Ua_2d_mat(:,:) = 0. 1438 polarity_2d_mat(:,:) = 0. 1439 1440 1441 END DO 1442 IF(lwp) WRITE(numout,*) "diaharm_fast: Finshed postprocessing 2d velocity tidal parameters" 1443 ENDIF 1444 1445 CALL FLUSH(numout) 1446 1447 IF (ln_ana_uv3d) THEN 1448 IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess baroclinic velocity tidal parameters" 1449 CALL FLUSH(numout) 1450 DO jh=1,nb_ana 1451 1452 1453 tmp_u_amp_3d_mat(:,:,:) = 0. 1454 tmp_v_amp_3d_mat(:,:,:) = 0. 1455 tmp_u_phi_3d_mat(:,:,:) = 0. 1456 tmp_v_phi_3d_mat(:,:,:) = 0. 1457 1458 a_u_3d_mat(:,:,:) = 0. 1459 b_u_3d_mat(:,:,:) = 0. 1460 a_v_3d_mat(:,:,:) = 0. 1461 b_v_3d_mat(:,:,:) = 0. 1462 1463 qmax_3d_mat(:,:,:) = 0. 1464 qmin_3d_mat(:,:,:) = 0. 1465 1466 ecc_3d_mat(:,:,:) = 0. 1467 thetamax_3d_mat(:,:,:) =0. 1468 thetamin_3d_mat(:,:,:) = 0. 1469 1470 Qc_3d_mat(:,:,:) = 0. 1471 Qac_3d_mat(:,:,:) = 0. 1472 gc_3d_mat(:,:,:) = 0. 1473 gac_3d_mat(:,:,:) = 0. 1474 1475 Phi_Ua_3d_mat(:,:,:) = 0. 1476 dir_Ua_3d_mat(:,:,:) = 0. 1477 polarity_3d_mat(:,:,:) = 0. 1478 1479 1480 DO jk=1,jpkm1 1481 !DO jj = 2, nlcj ! - 1 1482 ! DO ji = 2, nlci ! - 1 1483 1484 DO jj = 2, jpjm1 1485 DO ji = 2, jpim1 1486 1487 ! do jj=2,nlcj 1488 ! do ji=2,nlci 1489 !IF ((umask(ji,jj) + umask(ji-1,jj)) == 0 ) CYCLE 1490 !IF ((vmask(ji,jj) + vmask(ji,jj-1)) == 0 ) CYCLE 1491 1492 IF ( ((umask(ji,jj,jk) + umask(ji-1,jj,jk)) > 0 ) .AND. ((vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) > 0 ) ) THEN 1493 tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk)) 1494 tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) 1495 ! WORK ON THE WRAP AROUND 1496 !tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk)) 1497 !tmp_v_phi = ((phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) 1498 1499 if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then 1500 !tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk)) 1501 tmp_u_phi = atan2((sin(phi_u3d(jh,ji,jj,jk)) + sin(phi_u3d(jh,ji-1,jj,jk))),(cos(phi_u3d(jh,ji,jj,jk)) + cos(phi_u3d(jh,ji-1,jj,jk)))) 1502 else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then 1503 tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) 1504 else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then 1505 tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) 1506 else 1507 cycle 1508 end if 1509 1510 1511 if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then 1512 !tmp_v_phi = ((phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) 1513 tmp_v_phi = atan2((sin(phi_v3d(jh,ji,jj,jk)) + sin(phi_v3d(jh,ji,jj-1,jk))),(cos(phi_v3d(jh,ji,jj,jk)) + cos(phi_v3d(jh,ji,jj-1,jk)))) 1514 else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then 1515 tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) 1516 else if ( (vmask(ji,jj,jk) == 0) .AND. (vmask(ji,jj-1,jk) == 1)) then 1517 tmp_v_phi = (phi_v3d(ji,jj-1,jk)*vmask(ji,jj-1,jk)) 1518 else 1519 cycle 1520 end if 1521 1522 1523 ! do jj=1,nlcj 1524 ! do ji=1,nlci 1525 1526 ! tmp_u_amp = ((amp_u2d(jh,ji,jj)) + (amp_u2d(jh,ji-1,jj)))/(2.) 1527 ! tmp_v_amp = ((amp_v2d(jh,ji,jj)) + (amp_v2d(jh,ji,jj-1)))/(2.) 1528 ! ! WORK ON THE WRAP AROUND 1529 ! tmp_u_phi = ((phi_u2d(jh,ji,jj)) + (phi_u2d(jh,ji-1,jj)))/(2.) 1530 ! tmp_v_phi = ((phi_v2d(jh,ji,jj)) + (phi_v2d(jh,ji,jj-1)))/(2.) 1531 1532 1533 1534 ! tmp_u_amp = (amp_u2d(jh,ji,jj)) 1535 ! tmp_v_amp = (amp_v2d(jh,ji,jj)) 1536 ! ! WORK ON THE WRAP AROUND 1537 ! tmp_u_phi = (phi_u2d(jh,ji,jj)) 1538 ! tmp_v_phi = (phi_v2d(jh,ji,jj)) 1539 1540 1541 1542 a_u = tmp_U_amp * cos(tmp_U_phi) 1543 b_u = tmp_U_amp * sin(tmp_U_phi) 1544 a_v = tmp_V_amp * cos(tmp_V_phi) 1545 b_v = tmp_V_amp * sin(tmp_V_phi) 1546 1547 twodelta = atan2( (tmp_V_amp**2 * sin( 2*(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp**2 + tmp_V_amp**2 * cos( 2*(tmp_U_phi - tmp_V_phi) ) ) ) 1548 delta = twodelta/2. 1549 1550 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1551 1552 tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 1553 if (tmpreal < 0) tmpreal = 0 !CYCLE 1554 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) ) 1555 alpha2 = sqrt( tmpreal ) 1556 if (alpha2 < 0) alpha2 = 0 !CYCLE 1557 alpha= sqrt( alpha2 ) 1558 1559 1560 !major and minor axis of the ellipse 1561 qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) 1562 !tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1563 !qmin = 0 1564 !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1565 1566 tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 1567 if (tmpreal < 0) tmpreal = 0 !CYCLE 1568 !qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. 1569 qmin = sqrt( tmpreal ) ! but always positive. 1570 1571 !eccentricity of ellipse 1572 1573 tmpreal = (qmax + qmin) 1574 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE 1575 !ecc = (qmax - qmin)/(qmax + qmin) 1576 ecc = (qmax - qmin)/(tmpreal) 1577 ! Angle of major and minor ellipse 1578 thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) 1579 thetamin = thetamax + rpi/2. 1580 1581 1582 1583 ! Rotary current components: Pugh A3.10 1584 ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors 1585 ! so Qc = clockwise = anticyclonic = negative 1586 ! and Qac = anticlockwise = cyclonic = negative 1587 1588 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1589 if (tmpreal < 0) tmpreal = 0! CYCLE 1590 !Qc = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1591 Qc = 0.5*sqrt( tmpreal ) 1592 1593 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 1594 if (tmpreal < 0) tmpreal = 0! CYCLE 1595 !Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) ) 1596 Qac = 0.5*sqrt( tmpreal ) 1597 1598 1599 gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1600 gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1601 1602 !Pugh A3.2 1603 Phi_Ua = -0.5*(gac - gc) 1604 dir_Ua = 0.5*(gac + gc) ! positive from x axis 1605 1606 tmpreal = qmax 1607 !if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 1608 if (tmpreal == 0) then 1609 polarity = 0 1610 else 1611 polarity = (Qac - Qc)/qmax 1612 endif 1613 1614 1615 1616 tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp 1617 tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp 1618 tmp_u_phi_3d_mat(ji,jj,jk) = tmp_u_phi 1619 tmp_v_phi_3d_mat(ji,jj,jk) = tmp_v_phi 1620 1621 1622 a_u_3d_mat(ji,jj,jk) = a_u 1623 b_u_3d_mat(ji,jj,jk) = b_u 1624 a_v_3d_mat(ji,jj,jk) = a_v 1625 b_v_3d_mat(ji,jj,jk) = b_v 1626 1627 qmax_3d_mat(ji,jj,jk) = qmax 1628 qmin_3d_mat(ji,jj,jk) = qmin 1629 1630 ecc_3d_mat(ji,jj,jk) = ecc 1631 thetamax_3d_mat(ji,jj,jk) = thetamax 1632 thetamin_3d_mat(ji,jj,jk) = thetamin 1633 1634 Qc_3d_mat(ji,jj,jk) = Qc 1635 Qac_3d_mat(ji,jj,jk) = Qac 1636 gc_3d_mat(ji,jj,jk) = gc 1637 gac_3d_mat(ji,jj,jk) = gac 1638 1639 Phi_Ua_3d_mat(ji,jj,jk) = Phi_Ua 1640 dir_Ua_3d_mat(ji,jj,jk) = dir_Ua 1641 polarity_3d_mat(ji,jj,jk) = polarity 1642 1643 ENDIF 1644 END DO 1645 END DO 1646 END DO 1647 1648 1649 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d' 1650 IF( iom_use(TRIM(tmp_name)) ) THEN 1651 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1652 CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:)) 1653 ENDIF 1654 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d' 1655 IF( iom_use(TRIM(tmp_name)) ) THEN 1656 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1657 CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:)) 1658 ENDIF 1659 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d' 1660 IF( iom_use(TRIM(tmp_name)) ) THEN 1661 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1662 CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:)) 1663 ENDIF 1664 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d' 1665 IF( iom_use(TRIM(tmp_name)) ) THEN 1666 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1667 CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:)) 1668 ENDIF 1669 1670 1671 1672 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d' 1673 IF( iom_use(TRIM(tmp_name)) ) THEN 1674 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1675 CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:)) 1676 ENDIF 1677 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d' 1678 IF( iom_use(TRIM(tmp_name)) ) THEN 1679 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1680 CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:)) 1681 ENDIF 1682 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d' 1683 IF( iom_use(TRIM(tmp_name)) ) THEN 1684 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1685 CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:)) 1686 ENDIF 1687 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d' 1688 IF( iom_use(TRIM(tmp_name)) ) THEN 1689 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1690 CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:)) 1691 ENDIF 1692 1693 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d' 1694 IF( iom_use(TRIM(tmp_name)) ) THEN 1695 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1696 CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:)) 1697 ENDIF 1698 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d' 1699 IF( iom_use(TRIM(tmp_name)) ) THEN 1700 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1701 CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:)) 1702 ENDIF 1703 1704 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d' 1705 IF( iom_use(TRIM(tmp_name)) ) THEN 1706 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1707 CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:)) 1708 ENDIF 1709 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d' 1710 IF( iom_use(TRIM(tmp_name)) ) THEN 1711 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1712 CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:)) 1713 ENDIF 1714 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d' 1715 IF( iom_use(TRIM(tmp_name)) ) THEN 1716 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1717 CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:)) 1718 ENDIF 1719 1720 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d' 1721 IF( iom_use(TRIM(tmp_name)) ) THEN 1722 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1723 CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:)) 1724 ENDIF 1725 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d' 1726 IF( iom_use(TRIM(tmp_name)) ) THEN 1727 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1728 CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:)) 1729 ENDIF 1730 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d' 1731 IF( iom_use(TRIM(tmp_name)) ) THEN 1732 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1733 CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:)) 1734 ENDIF 1735 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d' 1736 IF( iom_use(TRIM(tmp_name)) ) THEN 1737 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1738 CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:)) 1739 ENDIF 1740 1741 1742 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d' 1743 IF( iom_use(TRIM(tmp_name)) ) THEN 1744 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1745 CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:)) 1746 ENDIF 1747 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d' 1748 IF( iom_use(TRIM(tmp_name)) ) THEN 1749 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1750 CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:)) 1751 ENDIF 1752 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d' 1753 IF( iom_use(TRIM(tmp_name)) ) THEN 1754 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1755 CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:)) 1756 ENDIF 1757 1758 tmp_u_amp_3d_mat(:,:,:) = 0. 1759 tmp_v_amp_3d_mat(:,:,:) = 0. 1760 tmp_u_phi_3d_mat(:,:,:) = 0. 1761 tmp_v_phi_3d_mat(:,:,:) = 0. 1762 1763 a_u_3d_mat(:,:,:) = 0. 1764 b_u_3d_mat(:,:,:) = 0. 1765 a_v_3d_mat(:,:,:) = 0. 1766 b_v_3d_mat(:,:,:) = 0. 1767 1768 qmax_3d_mat(:,:,:) = 0. 1769 qmin_3d_mat(:,:,:) = 0. 1770 1771 ecc_3d_mat(:,:,:) = 0. 1772 thetamax_3d_mat(:,:,:) =0. 1773 thetamin_3d_mat(:,:,:) = 0. 1774 1775 Qc_3d_mat(:,:,:) = 0. 1776 Qac_3d_mat(:,:,:) = 0. 1777 gc_3d_mat(:,:,:) = 0. 1778 gac_3d_mat(:,:,:) = 0. 1779 1780 Phi_Ua_3d_mat(:,:,:) = 0. 1781 dir_Ua_3d_mat(:,:,:) = 0. 1782 polarity_3d_mat(:,:,:) = 0. 1783 1784 1785 END DO 1786 1787 1788 1789 1790 1791 IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess 3d velocity tidal parameters" 1792 ENDIF 1793 1794 CALL FLUSH(numout) 1795 ENDIF 1796 1797 1798 CALL FLUSH(numout) 1244 1799 1245 1800 ! to output tidal parameters, u and v on t grid … … 1258 1813 1259 1814 1260 ! IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN 1261 1262 ! DEALLOCATE(amp_u2d, amp_v2d, phi_u2d, phi_v2d ) 1263 1264 1265 ! DEALLOCATE(tmp_u_amp_mat, tmp_v_amp_mat, tmp_u_phi_mat, tmp_v_phi_mat ) 1266 !! DEALLOCATE(a_u_mat, b_u_mat, a_v_mat, b_v_mat, qmax_mat, qmin_mat, ecc_mat ) 1267 !! DEALLOCATE(thetamax_mat, thetamin_mat, Qc_mat, Qac_mat, gc_mat, gac_mat ) 1268 !! DEALLOCATE(Phi_Ua_mat, dir_Ua_mat, polarity_mat ) 1269 1270 ! endif 1271 1272 ! IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 2d velocity tidal parameters" 1815 IF (ln_diaharm_postproc_vel) THEN 1816 IF (ln_ana_uvbar) THEN 1817 1818 DEALLOCATE(amp_u2d, amp_v2d, phi_u2d, phi_v2d ) 1819 1820 DEALLOCATE(tmp_u_amp_2d_mat, tmp_v_amp_2d_mat, tmp_u_phi_2d_mat, tmp_v_phi_2d_mat ) 1821 1822 DEALLOCATE(a_u_2d_mat, b_u_2d_mat, a_v_2d_mat, b_v_2d_mat ) 1823 DEALLOCATE(qmax_2d_mat, qmin_2d_mat, ecc_2d_mat ) 1824 DEALLOCATE(thetamax_2d_mat, thetamin_2d_mat, Qc_2d_mat, Qac_2d_mat) 1825 DEALLOCATE(gc_2d_mat, gac_2d_mat, Phi_Ua_2d_mat, dir_Ua_2d_mat) 1826 DEALLOCATE(polarity_2d_mat ) 1827 1828 ENDIF 1829 IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 2d velocity tidal parameters" 1830 1831 IF (ln_ana_uv3d) THEN 1832 1833 DEALLOCATE(amp_u3d, amp_v3d, phi_u3d, phi_v3d ) 1834 1835 DEALLOCATE(tmp_u_amp_3d_mat, tmp_v_amp_3d_mat, tmp_u_phi_3d_mat, tmp_v_phi_3d_mat ) 1836 1837 DEALLOCATE(a_u_3d_mat, b_u_3d_mat, a_v_3d_mat, b_v_3d_mat ) 1838 DEALLOCATE(qmax_3d_mat, qmin_3d_mat, ecc_3d_mat ) 1839 DEALLOCATE(thetamax_3d_mat, thetamin_3d_mat, Qc_3d_mat, Qac_3d_mat) 1840 DEALLOCATE(gc_3d_mat, gac_3d_mat, Phi_Ua_3d_mat, dir_Ua_3d_mat) 1841 DEALLOCATE(polarity_3d_mat ) 1842 1843 ENDIF 1844 IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 3d velocity tidal parameters" 1845 1846 ENDIF 1847 1273 1848 1274 1849 CALL FLUSH(numout)
Note: See TracChangeset
for help on using the changeset viewer.
