Changeset 15608 for NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src
- Timestamp:
- 2021-12-16T20:22:34+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
r15600 r15608 83 83 LOGICAL, PUBLIC :: ln_diaharm_fast 84 84 LOGICAL, PUBLIC :: ln_diaharm_postproc_vel 85 LOGICAL, PUBLIC :: ln_diaharm_verbose 85 86 86 87 … … 179 180 180 181 181 IF(lwp ) WRITE(numout,*) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)*86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp182 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)*86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp 182 183 183 184 IF( iom_use('tide_t') ) CALL iom_put( 'tide_t', sec2start ) … … 366 367 367 368 368 NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, & 369 & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily,ln_diaharm_postproc_vel 369 NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, & 370 & ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, & 371 & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily,ln_diaharm_postproc_vel, ln_diaharm_verbose 370 372 !!---------------------------------------------------------------------- 371 373 !JT … … 416 418 WRITE(numout,*) ' Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily 417 419 WRITE(numout,*) ' Number of Harmonics: nyear, nmonth = ', nyear, nmonth 418 WRITE(numout,*) ' Post-process velocity stats: ln_diaharm_postproc_vel = ', ln_diaharm_postproc_vel 420 WRITE(numout,*) ' Post-process velocity stats: ln_diaharm_postproc_vel = ',ln_diaharm_postproc_vel 421 WRITE(numout,*) ' Verbose: ln_diaharm_verbose = ', ln_diaharm_verbose 419 422 420 423 ENDIF … … 447 450 IF ( kt < 10 ) THEN 448 451 ln_diaharm_read_restart = .FALSE. 449 IF(lwp ) WRITE(numout,*) ' kt = ',kt450 IF(lwp ) WRITE(numout,*) ' kt < 10, so setting ln_diaharm_read_restart to .FALSE.'451 IF(lwp ) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart452 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' kt = ',kt 453 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' kt < 10, so setting ln_diaharm_read_restart to .FALSE.' 454 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart 452 455 ENDIF 453 456 … … 777 780 778 781 IF (ln_diaharm_postproc_vel) THEN 779 IF (ln_ana_uvbar) 782 IF (ln_ana_uvbar) THEN 780 783 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 784 783 785 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 786 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)) … … 791 793 792 794 793 IF (ln_ana_uv3d) 795 IF (ln_ana_uv3d) THEN 794 796 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 797 798 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 799 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)) … … 842 843 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 843 844 IF( iom_use(TRIM(tmp_name)) ) THEN 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"845 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D) 846 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr" 846 847 CALL iom_put( TRIM(tmp_name), h_out2D(:,:) ) 847 848 ELSE 848 IF(lwp ) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)849 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 849 850 ENDIF 850 851 851 852 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 852 853 IF( iom_use(TRIM(tmp_name)) ) THEN 853 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)854 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D) 854 855 CALL iom_put( TRIM(tmp_name), g_out2D(:,:) ) 855 856 ELSE 856 IF(lwp ) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)857 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 857 858 ENDIF 858 859 … … 902 903 tmp_name='TA_'//TRIM(suffix)//'_off' 903 904 IF( iom_use(TRIM(tmp_name)) ) THEN 904 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)905 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 905 906 CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid)) 906 907 ELSE 907 IF(lwp ) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)908 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 908 909 ENDIF 909 910 … … 953 954 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 954 955 IF( iom_use(TRIM(tmp_name)) ) THEN 955 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)956 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D) 956 957 CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) ) 957 958 ELSE 958 IF(lwp ) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)959 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 959 960 ENDIF 960 961 961 962 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 962 963 IF( iom_use(TRIM(tmp_name)) ) THEN 963 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)964 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D) 964 965 CALL iom_put(tmp_name, g_out3D(:,:,:) ) 965 966 ELSE 966 IF(lwp ) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)967 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 967 968 ENDIF 968 969 enddo ! jh970 971 suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )972 tmp_name='TA_'//TRIM(suffix)//'_off'973 IF( iom_use(TRIM(tmp_name)) ) THEN974 IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)975 CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))976 ELSE977 IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)978 ENDIF979 980 981 982 969 983 970 … … 1024 1011 CALL FLUSH(numout) 1025 1012 1013 1014 enddo ! jh 1015 1016 suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) ) 1017 tmp_name='TA_'//TRIM(suffix)//'_off' 1018 IF( iom_use(TRIM(tmp_name)) ) THEN 1019 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1020 CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid)) 1021 ELSE 1022 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 1023 ENDIF 1024 1026 1025 enddo ! jgrid 1027 1026 1028 1027 CALL FLUSH(numout) 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1028 1045 1029 … … 1078 1062 1079 1063 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 1064 DO jj = 1, nlcj !- 1 1065 DO ji = 1, nlci ! - 1 1087 1066 1088 1067 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 AROUND1092 !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 1068 1095 1069 if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then 1070 1071 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)) 1096 1072 !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 1073 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 1074 else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then 1075 tmp_u_amp = (amp_u2d(jh,ji,jj)*ssumask(ji,jj)) 1099 1076 tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj)) 1100 1077 else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then 1078 tmp_u_amp = (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) 1101 1079 tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) 1102 1080 else … … 1106 1084 1107 1085 if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then 1086 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)) 1108 1087 !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 1088 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 1089 else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then 1090 tmp_v_amp = (amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) 1111 1091 tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) 1112 1092 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)) 1093 tmp_v_amp = (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)) 1094 !tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)) 1095 tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)) 1114 1096 else 1115 1097 cycle 1116 1098 end if 1117 1099 1118 1119 ! do jj=1,nlcj1120 ! do ji=1,nlci1121 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 AROUND1125 ! 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 AROUND1133 ! 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) CYCLE1150 ! 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) CYCLE1152 ! alpha= sqrt( alpha2 )1153 1154 1155 ! !major and minor axis of the ellipse1156 ! 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)/21158 ! !qmin = 01159 ! !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)/21162 ! if (tmpreal < 0) CYCLE1163 ! qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive.1164 1165 ! !eccentricity of ellipse1166 ! tmpreal = (qmax + qmin)1167 ! if (tmpreal < 0) CYCLE1168 ! ecc = (qmax - qmin)/(qmax + qmin)1169 ! ! Angle of major and minor ellipse1170 ! 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.101176 ! ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors1177 ! ! so Qc = clockwise = anticyclonic = negative1178 ! ! and Qac = anticlockwise = cyclonic = negative1179 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) CYCLE1182 ! 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) CYCLE1186 ! 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.21193 ! Phi_Ua = -0.5*(gac - gc)1194 ! dir_Ua = 0.5*(gac + gc) ! positive from x axis1195 1196 ! tmpreal = qmax1197 ! if (tmpreal < 0) CYCLE1198 ! polarity = (Qac - Qc)/qmax1199 1100 1200 1101 a_u = tmp_U_amp * cos(tmp_U_phi) … … 1207 1108 1208 1109 !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 1110 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 1111 if (tmpreal < 0) tmpreal = 0 !CYCLE 1112 1212 1113 !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 1114 alpha2 = sqrt( tmpreal ) … … 1218 1119 !major and minor axis of the ellipse 1219 1120 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)/21221 !qmin = 01222 !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive.1223 1121 1224 1122 tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 … … 1228 1126 1229 1127 !eccentricity of ellipse 1230 1231 1128 tmpreal = (qmax + qmin) 1232 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE1129 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 1233 1130 !ecc = (qmax - qmin)/(qmax + qmin) 1234 1131 ecc = (qmax - qmin)/(tmpreal) 1132 1235 1133 ! Angle of major and minor ellipse 1236 1134 thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) … … 1245 1143 1246 1144 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 !CYCLE1145 if (tmpreal < 0) tmpreal = 0 !CYCLE 1248 1146 !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 1147 Qc = 0.5*sqrt( tmpreal ) 1250 1148 1251 1149 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 !CYCLE1150 if (tmpreal < 0) tmpreal = 0 !CYCLE 1253 1151 !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 1152 Qac = 0.5*sqrt( tmpreal ) 1255 1153 1256 1154 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 )) ) ) 1155 gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & 1156 & ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1157 gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & 1158 & ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1259 1159 1260 1160 !Pugh A3.2 … … 1270 1170 endif 1271 1171 1272 1273 1172 tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp 1274 1173 tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp … … 1305 1204 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar' 1306 1205 IF( iom_use(TRIM(tmp_name)) ) THEN 1307 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1206 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1308 1207 CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:)) 1309 1208 ENDIF 1310 1209 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar' 1311 1210 IF( iom_use(TRIM(tmp_name)) ) THEN 1312 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1211 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1313 1212 CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:)) 1314 1213 ENDIF 1315 1214 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar' 1316 1215 IF( iom_use(TRIM(tmp_name)) ) THEN 1317 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1216 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1318 1217 CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:)) 1319 1218 ENDIF 1320 1219 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar' 1321 1220 IF( iom_use(TRIM(tmp_name)) ) THEN 1322 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1221 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1323 1222 CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:)) 1324 1223 ENDIF … … 1328 1227 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar' 1329 1228 IF( iom_use(TRIM(tmp_name)) ) THEN 1330 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1229 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1331 1230 CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:)) 1332 1231 ENDIF 1333 1232 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar' 1334 1233 IF( iom_use(TRIM(tmp_name)) ) THEN 1335 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1234 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1336 1235 CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:)) 1337 1236 ENDIF 1338 1237 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar' 1339 1238 IF( iom_use(TRIM(tmp_name)) ) THEN 1340 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1239 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1341 1240 CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:)) 1342 1241 ENDIF 1343 1242 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar' 1344 1243 IF( iom_use(TRIM(tmp_name)) ) THEN 1345 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1244 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1346 1245 CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:)) 1347 1246 ENDIF … … 1349 1248 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar' 1350 1249 IF( iom_use(TRIM(tmp_name)) ) THEN 1351 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1250 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1352 1251 CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:)) 1353 1252 ENDIF 1354 1253 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar' 1355 1254 IF( iom_use(TRIM(tmp_name)) ) THEN 1356 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1255 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1357 1256 CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:)) 1358 1257 ENDIF … … 1360 1259 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar' 1361 1260 IF( iom_use(TRIM(tmp_name)) ) THEN 1362 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1261 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1363 1262 CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:)) 1364 1263 ENDIF 1365 1264 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar' 1366 1265 IF( iom_use(TRIM(tmp_name)) ) THEN 1367 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1266 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1368 1267 CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:)) 1369 1268 ENDIF 1370 1269 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar' 1371 1270 IF( iom_use(TRIM(tmp_name)) ) THEN 1372 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1271 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1373 1272 CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:)) 1374 1273 ENDIF … … 1376 1275 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar' 1377 1276 IF( iom_use(TRIM(tmp_name)) ) THEN 1378 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1277 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1379 1278 CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:)) 1380 1279 ENDIF 1381 1280 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar' 1382 1281 IF( iom_use(TRIM(tmp_name)) ) THEN 1383 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1282 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1384 1283 CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:)) 1385 1284 ENDIF 1386 1285 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar' 1387 1286 IF( iom_use(TRIM(tmp_name)) ) THEN 1388 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1287 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1389 1288 CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:)) 1390 1289 ENDIF 1391 1290 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar' 1392 1291 IF( iom_use(TRIM(tmp_name)) ) THEN 1393 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1292 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1394 1293 CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:)) 1395 1294 ENDIF … … 1398 1297 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar' 1399 1298 IF( iom_use(TRIM(tmp_name)) ) THEN 1400 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1299 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1401 1300 CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:)) 1402 1301 ENDIF 1403 1302 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar' 1404 1303 IF( iom_use(TRIM(tmp_name)) ) THEN 1405 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1304 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1406 1305 CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:)) 1407 1306 ENDIF 1408 1307 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar' 1409 1308 IF( iom_use(TRIM(tmp_name)) ) THEN 1410 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1309 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1411 1310 CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:)) 1412 1311 ENDIF … … 1478 1377 1479 1378 1480 1379 DO jk=1,jpkm1 1481 1380 !DO jj = 2, nlcj ! - 1 1482 1381 ! DO ji = 2, nlci ! - 1 1483 1382 1484 DO jj = 2, jpjm1 1485 DO ji = 2, jpim1 1383 ! DO jj = 1, jpjm1 #works 1384 ! DO ji = 1, jpim1 1385 1386 DO jj = 1, nlcj !- 1 1387 DO ji = 1, nlci ! - 1 1388 1486 1389 1487 1390 ! do jj=2,nlcj … … 1491 1394 1492 1395 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))1396 !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)) 1397 !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 1398 ! WORK ON THE WRAP AROUND 1496 1399 !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)) … … 1498 1401 1499 1402 if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then 1403 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)) 1500 1404 !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 1405 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 1406 else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then 1407 tmp_u_amp = (amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) 1503 1408 tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) 1504 1409 else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then 1410 tmp_u_amp = (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) 1505 1411 tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) 1506 1412 else 1507 1413 cycle 1508 end if 1414 end if 1509 1415 1510 1416 1511 1417 if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then 1418 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)) 1512 1419 !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 1420 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 1421 else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then 1422 tmp_v_amp = (amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) 1515 1423 tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) 1516 1424 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)) 1425 tmp_v_amp = (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)) 1426 tmp_v_phi = (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)) 1518 1427 else 1519 1428 cycle 1520 1429 end if 1521 1522 1523 ! do jj=1,nlcj1524 ! do ji=1,nlci1525 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 AROUND1529 ! 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 AROUND1537 ! tmp_u_phi = (phi_u2d(jh,ji,jj))1538 ! tmp_v_phi = (phi_v2d(jh,ji,jj))1539 1430 1540 1431 … … 1549 1440 1550 1441 !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 1442 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 1443 if (tmpreal < 0) tmpreal = 0 !CYCLE 1444 1554 1445 !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 1446 alpha2 = sqrt( tmpreal ) … … 1560 1451 !major and minor axis of the ellipse 1561 1452 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)/21563 !qmin = 01564 !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive.1565 1453 1566 1454 tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 … … 1570 1458 1571 1459 !eccentricity of ellipse 1572 1573 1460 tmpreal = (qmax + qmin) 1574 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE1461 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 1575 1462 !ecc = (qmax - qmin)/(qmax + qmin) 1576 1463 ecc = (qmax - qmin)/(tmpreal) 1464 1577 1465 ! Angle of major and minor ellipse 1578 1466 thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) … … 1587 1475 1588 1476 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 !CYCLE1477 if (tmpreal < 0) tmpreal = 0 !CYCLE 1590 1478 !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 1479 Qc = 0.5*sqrt( tmpreal ) 1592 1480 1593 1481 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 !CYCLE1482 if (tmpreal < 0) tmpreal = 0 !CYCLE 1595 1483 !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 1484 Qac = 0.5*sqrt( tmpreal ) 1597 1485 1598 1486 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 )) ) ) 1487 gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & 1488 & ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) 1489 gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & 1490 & ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) 1601 1491 1602 1492 !Pugh A3.2 … … 1613 1503 1614 1504 1615 1616 1505 tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp 1617 1506 tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp … … 1642 1531 1643 1532 ENDIF 1644 END DO 1645 END DO 1646 END DO 1533 END DO !ji 1534 END DO !jj 1535 END DO !jk 1647 1536 1648 1537 1649 1538 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d' 1650 1539 IF( iom_use(TRIM(tmp_name)) ) THEN 1651 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1540 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1652 1541 CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:)) 1653 1542 ENDIF 1654 1543 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d' 1655 1544 IF( iom_use(TRIM(tmp_name)) ) THEN 1656 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1545 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1657 1546 CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:)) 1658 1547 ENDIF 1659 1548 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d' 1660 1549 IF( iom_use(TRIM(tmp_name)) ) THEN 1661 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1550 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1662 1551 CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:)) 1663 1552 ENDIF 1664 1553 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d' 1665 1554 IF( iom_use(TRIM(tmp_name)) ) THEN 1666 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1555 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1667 1556 CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:)) 1668 1557 ENDIF … … 1672 1561 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d' 1673 1562 IF( iom_use(TRIM(tmp_name)) ) THEN 1674 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1563 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1675 1564 CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:)) 1676 1565 ENDIF 1677 1566 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d' 1678 1567 IF( iom_use(TRIM(tmp_name)) ) THEN 1679 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1568 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1680 1569 CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:)) 1681 1570 ENDIF 1682 1571 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d' 1683 1572 IF( iom_use(TRIM(tmp_name)) ) THEN 1684 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1573 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1685 1574 CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:)) 1686 1575 ENDIF 1687 1576 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d' 1688 1577 IF( iom_use(TRIM(tmp_name)) ) THEN 1689 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1578 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1690 1579 CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:)) 1691 1580 ENDIF … … 1693 1582 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d' 1694 1583 IF( iom_use(TRIM(tmp_name)) ) THEN 1695 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1584 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1696 1585 CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:)) 1697 1586 ENDIF 1698 1587 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d' 1699 1588 IF( iom_use(TRIM(tmp_name)) ) THEN 1700 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1589 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1701 1590 CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:)) 1702 1591 ENDIF … … 1704 1593 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d' 1705 1594 IF( iom_use(TRIM(tmp_name)) ) THEN 1706 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1595 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1707 1596 CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:)) 1708 1597 ENDIF 1709 1598 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d' 1710 1599 IF( iom_use(TRIM(tmp_name)) ) THEN 1711 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1600 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1712 1601 CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:)) 1713 1602 ENDIF 1714 1603 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d' 1715 1604 IF( iom_use(TRIM(tmp_name)) ) THEN 1716 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1605 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1717 1606 CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:)) 1718 1607 ENDIF … … 1720 1609 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d' 1721 1610 IF( iom_use(TRIM(tmp_name)) ) THEN 1722 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1611 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1723 1612 CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:)) 1724 1613 ENDIF 1725 1614 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d' 1726 1615 IF( iom_use(TRIM(tmp_name)) ) THEN 1727 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1616 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1728 1617 CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:)) 1729 1618 ENDIF 1730 1619 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d' 1731 1620 IF( iom_use(TRIM(tmp_name)) ) THEN 1732 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1621 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1733 1622 CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:)) 1734 1623 ENDIF 1735 1624 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d' 1736 1625 IF( iom_use(TRIM(tmp_name)) ) THEN 1737 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1626 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1738 1627 CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:)) 1739 1628 ENDIF … … 1742 1631 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d' 1743 1632 IF( iom_use(TRIM(tmp_name)) ) THEN 1744 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1633 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1745 1634 CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:)) 1746 1635 ENDIF 1747 1636 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d' 1748 1637 IF( iom_use(TRIM(tmp_name)) ) THEN 1749 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1638 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1750 1639 CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:)) 1751 1640 ENDIF 1752 1641 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d' 1753 1642 IF( iom_use(TRIM(tmp_name)) ) THEN 1754 IF(lwp ) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)1643 IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 1755 1644 CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:)) 1756 1645 ENDIF … … 1770 1659 1771 1660 ecc_3d_mat(:,:,:) = 0. 1772 thetamax_3d_mat(:,:,:) = 0.1661 thetamax_3d_mat(:,:,:) = 0. 1773 1662 thetamin_3d_mat(:,:,:) = 0. 1774 1663 … … 1783 1672 1784 1673 1785 END DO 1674 END DO !jh 1786 1675 1787 1676
Note: See TracChangeset
for help on using the changeset viewer.