New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15608 – NEMO

Changeset 15608


Ignore:
Timestamp:
2021-12-16T20:22:34+01:00 (2 years ago)
Author:
hadjt
Message:

DIA/diaharm_fast.F90 working

bug fixed diaharm_fast.F90

  • amp_u3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) was outside the harmonics loop, so only ran for M2
  • this also fixed the vertical stripes in the 3d output fields

Code also tidied up, and a namelist verbosity switch was added

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  
    8383   LOGICAL, PUBLIC :: ln_diaharm_fast 
    8484   LOGICAL, PUBLIC :: ln_diaharm_postproc_vel 
     85   LOGICAL, PUBLIC :: ln_diaharm_verbose 
    8586 
    8687 
    … …  
    179180 
    180181           
    181           IF(lwp) WRITE(numout,*) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)*86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp 
     182          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 
    182183 
    183184          IF( iom_use('tide_t') ) CALL iom_put( 'tide_t', sec2start ) 
    … …  
    366367       
    367368 
    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 
    370372      !!---------------------------------------------------------------------- 
    371373      !JT 
    … …  
    416418         WRITE(numout,*) '   Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily 
    417419         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 
    419422 
    420423      ENDIF 
    … …  
    447450      IF ( kt < 10 ) THEN 
    448451        ln_diaharm_read_restart = .FALSE. 
    449         IF(lwp) WRITE(numout,*) '   kt = ',kt 
    450         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_restart         
     452        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         
    452455      ENDIF 
    453456 
    … …  
    777780 
    778781      IF (ln_diaharm_postproc_vel)  THEN 
    779           IF (ln_ana_uvbar)  THEN 
     782          IF (ln_ana_uvbar) THEN 
    780783             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              
    783785             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)) 
    784786             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)) 
    … …  
    791793 
    792794 
    793           IF (ln_ana_uv3d)  THEN 
     795          IF (ln_ana_uv3d) THEN 
    794796             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              
    797798             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)) 
    798799             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)) 
    … …  
    842843             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 
    843844             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" 
    846847                CALL iom_put( TRIM(tmp_name), h_out2D(:,:) ) 
    847848             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) 
    849850             ENDIF 
    850851 
    851852             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 
    852853             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) 
    854855                CALL iom_put( TRIM(tmp_name), g_out2D(:,:) ) 
    855856             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) 
    857858             ENDIF 
    858859 
    … …  
    902903         tmp_name='TA_'//TRIM(suffix)//'_off' 
    903904         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) 
    905906            CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid)) 
    906907         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) 
    908909         ENDIF 
    909910 
    … …  
    953954             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) 
    954955             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) 
    956957                CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) ) 
    957958             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) 
    959960             ENDIF 
    960961 
    961962             tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) 
    962963             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) 
    964965                CALL iom_put(tmp_name, g_out3D(:,:,:) ) 
    965966             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) 
    967968             ENDIF 
    968  
    969           enddo             ! jh  
    970  
    971          suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) ) 
    972          tmp_name='TA_'//TRIM(suffix)//'_off' 
    973          IF( iom_use(TRIM(tmp_name)) )  THEN 
    974             IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) 
    975             CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid)) 
    976          ELSE 
    977             IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) 
    978          ENDIF 
    979  
    980  
    981  
    982969 
    983970 
    … …  
    10241011             CALL FLUSH(numout) 
    10251012 
     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 
    10261025      enddo                 ! jgrid 
    10271026 
    10281027     CALL FLUSH(numout) 
    1029  
    1030  
    1031  
    1032  
    1033  
    1034  
    1035  
    1036  
    1037  
    1038  
    1039  
    1040  
    1041  
    1042  
    1043  
    10441028 
    10451029 
    … …  
    10781062 
    10791063 
    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 
    10871066 
    10881067                        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)) 
    10941068 
    10951069                            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)) 
    10961072                              !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)) 
    10971073                              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)))) 
    10981074                            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)) 
    10991076                              tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj)) 
    11001077                            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)) 
    11011079                              tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) 
    11021080                            else  
    … …  
    11061084 
    11071085                            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)) 
    11081087                              !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)) 
    11091088                              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)))) 
    11101089                            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)) 
    11111091                              tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) 
    11121092                            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)) 
    11141096                            else  
    11151097                              cycle 
    11161098                            end if 
    11171099 
    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 
    11991100 
    12001101                            a_u = tmp_U_amp * cos(tmp_U_phi) 
    … …  
    12071108 
    12081109                            !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  
    12101110                            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))  
    12111111                            if (tmpreal < 0) tmpreal = 0 !CYCLE 
     1112 
    12121113                            !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))  ) 
    12131114                            alpha2 = sqrt( tmpreal ) 
    … …  
    12181119                            !major and minor axis of the ellipse 
    12191120                            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. 
    12231121 
    12241122                            tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 
    … …  
    12281126 
    12291127                            !eccentricity of ellipse 
    1230  
    12311128                            tmpreal =  (qmax + qmin) 
    1232                             if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE 
     1129                            if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 
    12331130                            !ecc = (qmax - qmin)/(qmax + qmin) 
    12341131                            ecc = (qmax - qmin)/(tmpreal) 
     1132 
    12351133                            ! Angle of major and minor ellipse 
    12361134                            thetamax = atan2((  tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta)   ) , ( tmp_U_amp * cos( delta) )  ) 
    … …  
    12451143 
    12461144                            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 
     1145                            if (tmpreal < 0) tmpreal = 0 !CYCLE 
    12481146                            !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))  ) 
    12491147                            Qc  = 0.5*sqrt( tmpreal ) 
    12501148 
    12511149                            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 
     1150                            if (tmpreal < 0) tmpreal = 0 !CYCLE 
    12531151                            !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))  ) 
    12541152                            Qac = 0.5*sqrt( tmpreal ) 
    12551153 
    12561154 
    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 ))  )  ) 
    12591159 
    12601160                            !Pugh A3.2 
    … …  
    12701170                            endif 
    12711171 
    1272  
    12731172                            tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp 
    12741173                            tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp 
    … …  
    13051204                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar' 
    13061205                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) 
    13081207                   CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:)) 
    13091208                ENDIF 
    13101209                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar' 
    13111210                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) 
    13131212                  CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:)) 
    13141213                ENDIF 
    13151214                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar' 
    13161215                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) 
    13181217                  CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:)) 
    13191218                ENDIF 
    13201219                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar' 
    13211220                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) 
    13231222                  CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:)) 
    13241223                ENDIF 
    … …  
    13281227                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar' 
    13291228                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) 
    13311230                   CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:)) 
    13321231                ENDIF 
    13331232                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar' 
    13341233                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) 
    13361235                  CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:)) 
    13371236                ENDIF 
    13381237                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar' 
    13391238                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) 
    13411240                  CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:)) 
    13421241                ENDIF 
    13431242                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar' 
    13441243                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) 
    13461245                  CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:)) 
    13471246                ENDIF 
    … …  
    13491248                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar' 
    13501249                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) 
    13521251                  CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:)) 
    13531252                ENDIF 
    13541253                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar' 
    13551254                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) 
    13571256                  CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:)) 
    13581257                ENDIF 
    … …  
    13601259                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar' 
    13611260                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) 
    13631262                  CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:)) 
    13641263                ENDIF 
    13651264                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar' 
    13661265                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) 
    13681267                  CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:)) 
    13691268                ENDIF 
    13701269                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar' 
    13711270                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) 
    13731272                  CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:)) 
    13741273                ENDIF 
    … …  
    13761275                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar' 
    13771276                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) 
    13791278                  CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:)) 
    13801279                ENDIF 
    13811280                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar' 
    13821281                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) 
    13841283                  CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:)) 
    13851284                ENDIF 
    13861285                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar' 
    13871286                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) 
    13891288                  CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:)) 
    13901289                ENDIF 
    13911290                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar' 
    13921291                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) 
    13941293                  CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:)) 
    13951294                ENDIF 
    … …  
    13981297                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar' 
    13991298                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) 
    14011300                  CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:)) 
    14021301                ENDIF 
    14031302                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar' 
    14041303                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) 
    14061305                  CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:)) 
    14071306                ENDIF 
    14081307                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar' 
    14091308                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) 
    14111310                  CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:)) 
    14121311                ENDIF 
    … …  
    14781377 
    14791378 
    1480                   DO jk=1,jpkm1 
     1379                DO jk=1,jpkm1 
    14811380                     !DO jj = 2, nlcj ! - 1 
    14821381                     !   DO ji = 2, nlci ! - 1 
    14831382 
    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 
    14861389 
    14871390        !             do jj=2,nlcj 
    … …  
    14911394 
    14921395                            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)) 
    14951398                                ! WORK ON THE WRAP AROUND 
    14961399                                !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)) 
    … …  
    14981401 
    14991402                                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)) 
    15001404                                  !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)) 
    15011405                                  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)))) 
    15021406                                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))  
    15031408                                  tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) 
    15041409                                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)) 
    15051411                                  tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) 
    15061412                                else  
    15071413                                  cycle 
    1508                                 end if 
     1414                                end if   
    15091415 
    15101416 
    15111417                                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)) 
    15121419                                  !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)) 
    15131420                                  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)))) 
    15141421                                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)) 
    15151423                                  tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) 
    15161424                                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)) 
    15181427                                else  
    15191428                                  cycle 
    15201429                                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))  
    15391430 
    15401431 
    … …  
    15491440 
    15501441                                !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  
    15521442                                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))  
    15531443                                if (tmpreal < 0) tmpreal = 0 !CYCLE 
     1444 
    15541445                                !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))  ) 
    15551446                                alpha2 = sqrt( tmpreal ) 
    … …  
    15601451                                !major and minor axis of the ellipse 
    15611452                                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. 
    15651453 
    15661454                                tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 
    … …  
    15701458 
    15711459                                !eccentricity of ellipse 
    1572  
    15731460                                tmpreal =  (qmax + qmin) 
    1574                                 if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE 
     1461                                if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE 
    15751462                                !ecc = (qmax - qmin)/(qmax + qmin) 
    15761463                                ecc = (qmax - qmin)/(tmpreal) 
     1464 
    15771465                                ! Angle of major and minor ellipse 
    15781466                                thetamax = atan2((  tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta)   ) , ( tmp_U_amp * cos( delta) )  ) 
    … …  
    15871475 
    15881476                                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 
     1477                                if (tmpreal < 0) tmpreal = 0 !CYCLE 
    15901478                                !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))  ) 
    15911479                                Qc  = 0.5*sqrt( tmpreal ) 
    15921480 
    15931481                                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 
     1482                                if (tmpreal < 0) tmpreal = 0 !CYCLE 
    15951483                                !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))  ) 
    15961484                                Qac = 0.5*sqrt( tmpreal ) 
    15971485 
    15981486 
    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 ))  )  ) 
    16011491 
    16021492                                !Pugh A3.2 
    … …  
    16131503 
    16141504 
    1615  
    16161505                                tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp 
    16171506                                tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp 
    … …  
    16421531 
    16431532                            ENDIF 
    1644                         END DO 
    1645                      END DO 
    1646                  END DO 
     1533                        END DO !ji 
     1534                     END DO    !jj 
     1535                 END DO        !jk 
    16471536 
    16481537 
    16491538                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d' 
    16501539                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) 
    16521541                   CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:)) 
    16531542                ENDIF 
    16541543                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d' 
    16551544                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) 
    16571546                  CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:)) 
    16581547                ENDIF 
    16591548                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d' 
    16601549                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) 
    16621551                  CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:)) 
    16631552                ENDIF 
    16641553                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d' 
    16651554                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) 
    16671556                  CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:)) 
    16681557                ENDIF 
    … …  
    16721561                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d' 
    16731562                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) 
    16751564                   CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:)) 
    16761565                ENDIF 
    16771566                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d' 
    16781567                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) 
    16801569                  CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:)) 
    16811570                ENDIF 
    16821571                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d' 
    16831572                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) 
    16851574                  CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:)) 
    16861575                ENDIF 
    16871576                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d' 
    16881577                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) 
    16901579                  CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:)) 
    16911580                ENDIF 
    … …  
    16931582                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d' 
    16941583                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) 
    16961585                  CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:)) 
    16971586                ENDIF 
    16981587                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d' 
    16991588                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) 
    17011590                  CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:)) 
    17021591                ENDIF 
    … …  
    17041593                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d' 
    17051594                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) 
    17071596                  CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:)) 
    17081597                ENDIF 
    17091598                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d' 
    17101599                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) 
    17121601                  CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:)) 
    17131602                ENDIF 
    17141603                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d' 
    17151604                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) 
    17171606                  CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:)) 
    17181607                ENDIF 
    … …  
    17201609                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d' 
    17211610                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) 
    17231612                  CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:)) 
    17241613                ENDIF 
    17251614                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d' 
    17261615                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) 
    17281617                  CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:)) 
    17291618                ENDIF 
    17301619                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d' 
    17311620                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) 
    17331622                  CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:)) 
    17341623                ENDIF 
    17351624                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d' 
    17361625                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) 
    17381627                  CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:)) 
    17391628                ENDIF 
    … …  
    17421631                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d' 
    17431632                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) 
    17451634                  CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:)) 
    17461635                ENDIF 
    17471636                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d' 
    17481637                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) 
    17501639                  CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:)) 
    17511640                ENDIF 
    17521641                tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d' 
    17531642                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) 
    17551644                  CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:)) 
    17561645                ENDIF 
    … …  
    17701659 
    17711660                ecc_3d_mat(:,:,:) = 0. 
    1772                 thetamax_3d_mat(:,:,:) =0. 
     1661                thetamax_3d_mat(:,:,:) = 0. 
    17731662                thetamin_3d_mat(:,:,:) = 0. 
    17741663 
    … …  
    17831672 
    17841673 
    1785              END DO 
     1674             END DO    !jh 
    17861675 
    17871676 
Note: See TracChangeset for help on using the changeset viewer.