| 1 | MODULE diaharm_fast |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE example *** |
|---|
| 4 | !! Ocean physics: On line harmonic analyser |
|---|
| 5 | !! |
|---|
| 6 | !!===================================================================== |
|---|
| 7 | |
|---|
| 8 | !!---------------------------------------------------------------------- |
|---|
| 9 | !! : Calculate harmonic analysis |
|---|
| 10 | !!---------------------------------------------------------------------- |
|---|
| 11 | !! harm_ana : |
|---|
| 12 | !! harm_ana_init : |
|---|
| 13 | !! NB: 2017-12 : add 3D harmonic analysis of velocities |
|---|
| 14 | !! integration of Maria Luneva's development |
|---|
| 15 | !! 'key_3Ddiaharm' |
|---|
| 16 | !!---------------------------------------------------------------------- |
|---|
| 17 | |
|---|
| 18 | USE oce ! ocean dynamics and tracers |
|---|
| 19 | USE dom_oce ! ocean space and time domain |
|---|
| 20 | USE iom |
|---|
| 21 | USE in_out_manager ! I/O units |
|---|
| 22 | USE phycst ! physical constants |
|---|
| 23 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
|---|
| 24 | USE bdy_oce ! ocean open boundary conditions |
|---|
| 25 | USE bdytides ! tidal bdy forcing |
|---|
| 26 | USE daymod ! calendar |
|---|
| 27 | USE tideini |
|---|
| 28 | USE restart |
|---|
| 29 | USE ioipsl, ONLY : ju2ymds ! for calendar |
|---|
| 30 | ! |
|---|
| 31 | ! |
|---|
| 32 | USE timing ! preformance summary |
|---|
| 33 | USE zdf_oce |
|---|
| 34 | |
|---|
| 35 | IMPLICIT NONE |
|---|
| 36 | PRIVATE |
|---|
| 37 | |
|---|
| 38 | !! * Routine accessibility |
|---|
| 39 | PUBLIC dia_harm_fast ! routine called in step.F90 module |
|---|
| 40 | LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .TRUE. ! to be run or not |
|---|
| 41 | LOGICAL, PUBLIC :: lk_diaharm_2D ! = .TRUE. ! to run 2d |
|---|
| 42 | LOGICAL, PUBLIC :: lk_diaharm_3D ! = .TRUE. ! to run 3d |
|---|
| 43 | |
|---|
| 44 | !! * Module variables |
|---|
| 45 | INTEGER, PARAMETER :: nharm_max = jpmax_harmo ! max number of harmonics to be analysed |
|---|
| 46 | INTEGER, PARAMETER :: nhm_max = 2*nharm_max+1 |
|---|
| 47 | INTEGER, PARAMETER :: nvab = 2 ! number of 3D variables |
|---|
| 48 | INTEGER :: nharm |
|---|
| 49 | INTEGER :: nhm |
|---|
| 50 | INTEGER :: & !!! ** toto namelist (namtoto) ** |
|---|
| 51 | nflag = 1 ! default value of nflag |
|---|
| 52 | REAL(wp), DIMENSION(nharm_max) :: & |
|---|
| 53 | om_tide ! tidal frequencies ( rads/sec) |
|---|
| 54 | REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:) :: & |
|---|
| 55 | bzz,c,x ! work arrays |
|---|
| 56 | REAL(wp) :: cca,ssa,zm,bt,dd_cumul |
|---|
| 57 | ! |
|---|
| 58 | REAL(wp), PUBLIC :: fjulday_startharm !: Julian Day since start of harmonic analysis |
|---|
| 59 | REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana |
|---|
| 60 | REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a |
|---|
| 61 | ! |
|---|
| 62 | INTEGER :: nvar_2d, nvar_3d !: number of 2d and 3d variables to analyse |
|---|
| 63 | INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d |
|---|
| 64 | |
|---|
| 65 | ! Name of variables used in the restart |
|---|
| 66 | CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/) |
|---|
| 67 | CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/) |
|---|
| 68 | ! |
|---|
| 69 | REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,: ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D |
|---|
| 70 | REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D |
|---|
| 71 | ! |
|---|
| 72 | REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) :: g_out2D,h_out2D ! arrays for output |
|---|
| 73 | REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: g_out3D,h_out3D ! arrays for 3D output |
|---|
| 74 | ! |
|---|
| 75 | ! NAMELIST |
|---|
| 76 | LOGICAL, PUBLIC :: ln_diaharm_store !: =T Stores data for harmonic Analysis |
|---|
| 77 | LOGICAL, PUBLIC :: ln_diaharm_compute !: =T Compute harmonic Analysis |
|---|
| 78 | LOGICAL, PUBLIC :: ln_diaharm_read_restart !: =T Read restart from a previous run |
|---|
| 79 | !JT |
|---|
| 80 | LOGICAL, PUBLIC :: ln_diaharm_multiyear !: =T Read restart from a previous run |
|---|
| 81 | INTEGER, PUBLIC :: nn_diaharm_multiyear !: =T Read restart from a previous run |
|---|
| 82 | LOGICAL, PUBLIC :: ln_diaharm_update_nodal_daily !: =T update the nodes every day |
|---|
| 83 | LOGICAL, PUBLIC :: ln_diaharm_fast |
|---|
| 84 | LOGICAL, PUBLIC :: ln_diaharm_postproc_vel |
|---|
| 85 | LOGICAL, PUBLIC :: ln_diaharm_verbose |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | !JT |
|---|
| 89 | LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d |
|---|
| 90 | INTEGER :: nb_ana ! Number of harmonics to analyse |
|---|
| 91 | CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...) |
|---|
| 92 | INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_all ! INDEX within the full set of constituents (tide.h90) |
|---|
| 93 | INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_sub ! INDEX within the subset of constituents pass in input |
|---|
| 94 | |
|---|
| 95 | !! * Substitutions |
|---|
| 96 | |
|---|
| 97 | !!---------------------------------------------------------------------- |
|---|
| 98 | !! OPA 9.0 , LOCEAN-IPSL (2005) |
|---|
| 99 | !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005) |
|---|
| 100 | !! or TOP 1.0 , LOCEAN-IPSL (2005) |
|---|
| 101 | !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $ |
|---|
| 102 | !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt |
|---|
| 103 | !!---------------------------------------------------------------------- |
|---|
| 104 | |
|---|
| 105 | CONTAINS |
|---|
| 106 | |
|---|
| 107 | SUBROUTINE dia_harm_fast( kt ) |
|---|
| 108 | !!---------------------------------------------------------------------- |
|---|
| 109 | !! *** ROUTINE harm_ana *** |
|---|
| 110 | !! |
|---|
| 111 | !! ** Purpose : Harmonic analyser |
|---|
| 112 | !! |
|---|
| 113 | !! ** Method : |
|---|
| 114 | !! |
|---|
| 115 | !! ** Action : - first action (share memory array/varible modified |
|---|
| 116 | !! in this routine |
|---|
| 117 | !! - second action ..... |
|---|
| 118 | !! - ..... |
|---|
| 119 | !! |
|---|
| 120 | !! References : |
|---|
| 121 | !! Give references if exist otherwise suppress these lines |
|---|
| 122 | !! |
|---|
| 123 | !! History : |
|---|
| 124 | !! 9.0 ! 03-08 (Autor Names) Original code |
|---|
| 125 | !! ! 02-08 (Author names) brief description of modifications |
|---|
| 126 | !!---------------------------------------------------------------------- |
|---|
| 127 | !! * Modules used |
|---|
| 128 | |
|---|
| 129 | !! * arguments |
|---|
| 130 | INTEGER, INTENT( in ) :: & |
|---|
| 131 | kt ! describe it!!! |
|---|
| 132 | |
|---|
| 133 | !! * local declarations |
|---|
| 134 | INTEGER :: ji, jk, jj ! dummy loop arguments |
|---|
| 135 | INTEGER :: jh, i1, i2, jgrid |
|---|
| 136 | INTEGER :: j2d, j3d |
|---|
| 137 | REAL(WP) :: sec2start,sec2start_old |
|---|
| 138 | CHARACTER (len=40) :: tmp_name |
|---|
| 139 | !!-------------------------------------------------------------------- |
|---|
| 140 | |
|---|
| 141 | !JT IF( nn_timing == 1 ) CALL timing_start( 'dia_harm_fast' ) |
|---|
| 142 | IF( ln_timing ) CALL timing_start( 'dia_harm_fast' ) |
|---|
| 143 | IF( kt == nit000 ) CALL harm_ana_init(kt) ! Initialization (first time-step only) |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | IF ( ln_diaharm_update_nodal_daily ) THEN |
|---|
| 147 | !IF (MOD(kt,nint(86400./rdt)) == 0) THEN |
|---|
| 148 | |
|---|
| 149 | |
|---|
| 150 | IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN |
|---|
| 151 | DO jh = 1, nb_ana |
|---|
| 152 | !JT anau(jh) = 3.141579*utide ( ntide_sub(jh) )/180. |
|---|
| 153 | !JT anav(jh) = 3.141579*v0tide( ntide_sub(jh) )/180. |
|---|
| 154 | anau(jh) = utide ( ntide_sub(jh) ) |
|---|
| 155 | anav(jh) = v0tide( ntide_sub(jh) ) |
|---|
| 156 | anaf(jh) = ftide ( ntide_sub(jh) ) |
|---|
| 157 | END DO |
|---|
| 158 | |
|---|
| 159 | IF(lwp) WRITE(numout,*) |
|---|
| 160 | IF(lwp) WRITE(numout,*) 'harm_ana : update nodes?',ln_diaharm_update_nodal_daily |
|---|
| 161 | IF(lwp) WRITE(numout,*) 'harm_ana : date, time',ndastp, nhour, nminute |
|---|
| 162 | |
|---|
| 163 | ENDIF |
|---|
| 164 | ENDIF |
|---|
| 165 | |
|---|
| 166 | !IF( iom_use('diaharm_fast_eg_ssh_ts') ) THEN |
|---|
| 167 | ! IF (lwp) CALL iom_put( 'diaharm_fast_eg_ssh_ts', sshn(5,5) ) |
|---|
| 168 | !ENDIF |
|---|
| 169 | |
|---|
| 170 | IF ( ln_diaharm_fast .and. ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN |
|---|
| 171 | |
|---|
| 172 | ! this bit done every time step |
|---|
| 173 | nhm=2*nb_ana+1 |
|---|
| 174 | c(1) = 1.0 |
|---|
| 175 | |
|---|
| 176 | sec2start_old = nint( (fjulday-fjulday_startharm)*86400._wp ) |
|---|
| 177 | sec2start = nsec_day - NINT(0.5_wp * rdt) |
|---|
| 178 | |
|---|
| 179 | sec2start = adatrj * 86400._wp |
|---|
| 180 | |
|---|
| 181 | |
|---|
| 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 |
|---|
| 183 | |
|---|
| 184 | IF( iom_use('tide_t') ) CALL iom_put( 'tide_t', sec2start ) |
|---|
| 185 | |
|---|
| 186 | |
|---|
| 187 | |
|---|
| 188 | |
|---|
| 189 | |
|---|
| 190 | !IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm |
|---|
| 191 | |
|---|
| 192 | DO jh=1,nb_ana |
|---|
| 193 | c(2*jh ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) ) |
|---|
| 194 | c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) ) |
|---|
| 195 | |
|---|
| 196 | c(2*jh ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) ) |
|---|
| 197 | c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) ) |
|---|
| 198 | |
|---|
| 199 | |
|---|
| 200 | |
|---|
| 201 | !IF(lwp) WRITE(numout,*) 'diaharm_fast: analwave,',kt,tname(jh),sec2start,sec2start/3600.,sec2start_old,sec2start_old/3600,c(2*jh),c(2*jh+1),om_tide(jh),anau(jh),anav(jh) |
|---|
| 202 | ENDDO |
|---|
| 203 | |
|---|
| 204 | !IF(lwp) WRITE(numout,*) "c init", c, "c end", sec2start, om_tide(1), anau(1), anav(1),"end nodal" |
|---|
| 205 | |
|---|
| 206 | |
|---|
| 207 | ! CUMULATE |
|---|
| 208 | DO ji=1,jpi ! loop lon |
|---|
| 209 | DO jj=1,jpj ! loop lat |
|---|
| 210 | DO jh=1,nhm ! loop harmonic |
|---|
| 211 | |
|---|
| 212 | DO j2d=1,nvar_2d |
|---|
| 213 | IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj) ! analysis elevation |
|---|
| 214 | IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj) ! analysis depth average velocities |
|---|
| 215 | IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj) |
|---|
| 216 | !JT IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction |
|---|
| 217 | !JT IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj) |
|---|
| 218 | g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul |
|---|
| 219 | ENDDO |
|---|
| 220 | |
|---|
| 221 | DO j3d=1,nvar_3d |
|---|
| 222 | DO jk=1,jpkm1 |
|---|
| 223 | IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) * rhd(ji,jj,jk) * tmask(ji,jj,jk) |
|---|
| 224 | IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk) |
|---|
| 225 | IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk) |
|---|
| 226 | IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) * wn(ji,jj,jk) * wmask(ji,jj,jk) |
|---|
| 227 | g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul |
|---|
| 228 | ENDDO |
|---|
| 229 | ENDDO |
|---|
| 230 | |
|---|
| 231 | ENDDO ! end loop harmonic |
|---|
| 232 | ENDDO ! end loop lat |
|---|
| 233 | ENDDO ! end loop lon |
|---|
| 234 | |
|---|
| 235 | ! Compute nodal factor cumulative cross-product |
|---|
| 236 | DO i1=1,nhm |
|---|
| 237 | DO i2=1,nhm |
|---|
| 238 | cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2) |
|---|
| 239 | ENDDO |
|---|
| 240 | ENDDO |
|---|
| 241 | |
|---|
| 242 | ! Output RESTART |
|---|
| 243 | IF( kt == nitrst ) THEN |
|---|
| 244 | CALL harm_rst_write(kt) ! Dump out data for a restarted run |
|---|
| 245 | ENDIF |
|---|
| 246 | |
|---|
| 247 | ! At End of run |
|---|
| 248 | IF ( kt == nitend ) THEN |
|---|
| 249 | |
|---|
| 250 | IF(lwp) WRITE(numout,*) |
|---|
| 251 | IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run' |
|---|
| 252 | IF(lwp) WRITE(numout,*) '~~~~~~~~~' |
|---|
| 253 | |
|---|
| 254 | IF( ln_diaharm_compute ) THEN |
|---|
| 255 | |
|---|
| 256 | ! INITIALISE TABLE TO 0 |
|---|
| 257 | IF ( nvar_2d .gt. 0 ) THEN |
|---|
| 258 | g_cosamp2D = 0.0_wp |
|---|
| 259 | g_sinamp2D = 0.0_wp |
|---|
| 260 | ENDIF |
|---|
| 261 | IF ( nvar_3d .gt. 0 ) THEN |
|---|
| 262 | g_cosamp3D = 0.0_wp |
|---|
| 263 | g_sinamp3D = 0.0_wp |
|---|
| 264 | ENDIF |
|---|
| 265 | |
|---|
| 266 | ! FIRST OUTPUT 2D VARIABLES |
|---|
| 267 | DO jgrid=1,nvar_2d ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically |
|---|
| 268 | DO ji=1,jpi ! loop lon |
|---|
| 269 | DO jj=1,jpj ! loop lat |
|---|
| 270 | bt = 1.0_wp; bzz(:) = 0.0_wp |
|---|
| 271 | DO jh=1,nhm ! loop harmonic |
|---|
| 272 | bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid) |
|---|
| 273 | bt = bt*bzz(jh) |
|---|
| 274 | ENDDO |
|---|
| 275 | ! Copy back original cumulated nodal factor |
|---|
| 276 | a(:,:) = cc(:,:) |
|---|
| 277 | ! now do gaussian elimination of the system |
|---|
| 278 | ! a * x = b |
|---|
| 279 | ! the matrix x is (a0,a1,b1,a2,b2 ...) |
|---|
| 280 | ! the matrix a and rhs b solved here for x |
|---|
| 281 | x=0.0_wp |
|---|
| 282 | IF(bt.ne.0.) THEN |
|---|
| 283 | CALL gelim( a, bzz, x, nhm ) |
|---|
| 284 | ! Backup output in variables |
|---|
| 285 | DO jh=1,nb_ana |
|---|
| 286 | g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2 ) |
|---|
| 287 | g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1) |
|---|
| 288 | ENDDO |
|---|
| 289 | g_cosamp2D ( 0,ji,jj,jgrid) = x(1) |
|---|
| 290 | g_sinamp2D ( 0,ji,jj,jgrid) = 0.0_wp |
|---|
| 291 | ENDIF ! bt.ne.0. |
|---|
| 292 | ENDDO ! jj |
|---|
| 293 | ENDDO ! ji |
|---|
| 294 | ENDDO ! jgrid |
|---|
| 295 | |
|---|
| 296 | ! SECOND OUTPUT 3D VARIABLES |
|---|
| 297 | DO jgrid=1,nvar_3d ! loop number of 3d variables rho, U, V, W |
|---|
| 298 | DO jk=1,jpkm1 ! loop over vertical level |
|---|
| 299 | DO ji=1,jpi ! loop over lon |
|---|
| 300 | DO jj=1,jpj ! loop over lat |
|---|
| 301 | bt = 1.0_wp; bzz(:) = 0.0_wp |
|---|
| 302 | DO jh=1,nhm |
|---|
| 303 | bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid) |
|---|
| 304 | bt = bt*bzz(jh) |
|---|
| 305 | ENDDO |
|---|
| 306 | ! Copy back original cumulated nodal factor |
|---|
| 307 | a(:,:) = cc(:,:) |
|---|
| 308 | ! now do gaussian elimination of the system |
|---|
| 309 | ! a * x = b |
|---|
| 310 | ! the matrix x is (a0,a1,b1,a2,b2 ...) |
|---|
| 311 | ! the matrix a and rhs b solved here for x |
|---|
| 312 | x=0.0_wp |
|---|
| 313 | IF(bt.ne.0.) THEN |
|---|
| 314 | CALL gelim( a, bzz, x, nhm ) |
|---|
| 315 | ! Backup output in variables |
|---|
| 316 | DO jh=1,nb_ana |
|---|
| 317 | g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2 ) |
|---|
| 318 | g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1) |
|---|
| 319 | ENDDO |
|---|
| 320 | g_cosamp3D ( 0,ji,jj,jk,jgrid) = x(1) |
|---|
| 321 | g_sinamp3D ( 0,ji,jj,jk,jgrid) = 0.0_wp |
|---|
| 322 | ENDIF ! bt.ne.0. |
|---|
| 323 | ENDDO ! jj |
|---|
| 324 | ENDDO ! ji |
|---|
| 325 | ENDDO ! jk |
|---|
| 326 | ENDDO ! jgrid |
|---|
| 327 | |
|---|
| 328 | CALL harm_ana_out ! output analysis (last time step) |
|---|
| 329 | |
|---|
| 330 | ELSE ! ln_harmana_compute = False |
|---|
| 331 | IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step" |
|---|
| 332 | |
|---|
| 333 | ENDIF ! ln_harmana_compute |
|---|
| 334 | ENDIF ! kt == nitend |
|---|
| 335 | |
|---|
| 336 | ENDIF |
|---|
| 337 | |
|---|
| 338 | !JT IF( nn_timing == 1 ) CALL timing_stop( 'dia_harm_fast' ) |
|---|
| 339 | IF( ln_timing ) CALL timing_stop( 'dia_harm_fast' ) |
|---|
| 340 | |
|---|
| 341 | END SUBROUTINE dia_harm_fast |
|---|
| 342 | |
|---|
| 343 | SUBROUTINE harm_ana_init( kt ) !JT |
|---|
| 344 | !!---------------------------------------------------------------------- |
|---|
| 345 | !! *** ROUTINE harm_ana_init *** |
|---|
| 346 | !! |
|---|
| 347 | !! ** Purpose : initialization of .... |
|---|
| 348 | !! |
|---|
| 349 | !! ** Method : blah blah blah ... |
|---|
| 350 | !! |
|---|
| 351 | !! ** input : Namlist namexa |
|---|
| 352 | !! |
|---|
| 353 | !! ** Action : ... |
|---|
| 354 | !! |
|---|
| 355 | !! history : |
|---|
| 356 | !! 9.0 ! 03-08 (Autor Names) Original code |
|---|
| 357 | !!---------------------------------------------------------------------- |
|---|
| 358 | !! * local declarations |
|---|
| 359 | |
|---|
| 360 | INTEGER, INTENT(in) :: kt ! ocean time-step !JT |
|---|
| 361 | !! |
|---|
| 362 | INTEGER :: ji, jk, jh ! dummy loop indices |
|---|
| 363 | INTEGER :: ios ! Local integer output status for namelist read |
|---|
| 364 | INTEGER :: k2d, k3d ! dummy number of analysis |
|---|
| 365 | |
|---|
| 366 | |
|---|
| 367 | |
|---|
| 368 | |
|---|
| 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 |
|---|
| 372 | !!---------------------------------------------------------------------- |
|---|
| 373 | !JT |
|---|
| 374 | ln_diaharm_fast = .FALSE. |
|---|
| 375 | ln_diaharm_multiyear = .FALSE. |
|---|
| 376 | nn_diaharm_multiyear = 20 |
|---|
| 377 | !JT |
|---|
| 378 | lk_diaharm_2D = .TRUE. ! to run 2d |
|---|
| 379 | lk_diaharm_3D = .TRUE. ! to run 3d |
|---|
| 380 | |
|---|
| 381 | ln_diaharm_store = .TRUE. |
|---|
| 382 | |
|---|
| 383 | ln_diaharm_postproc_vel = .FALSE. |
|---|
| 384 | |
|---|
| 385 | IF(lwp) WRITE(numout,*) |
|---|
| 386 | IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides' |
|---|
| 387 | IF(lwp) WRITE(numout,*) '~~~~~~~~~' |
|---|
| 388 | |
|---|
| 389 | ! GET NAMELIST DETAILS |
|---|
| 390 | REWIND( numnam_ref ) ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis |
|---|
| 391 | READ ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901) |
|---|
| 392 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' ) |
|---|
| 393 | |
|---|
| 394 | REWIND( numnam_cfg ) ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis |
|---|
| 395 | READ ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 ) |
|---|
| 396 | 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' ) |
|---|
| 397 | IF(lwm) WRITE ( numond, nam_diaharm_fast ) |
|---|
| 398 | |
|---|
| 399 | |
|---|
| 400 | ! |
|---|
| 401 | IF(lwp) THEN |
|---|
| 402 | WRITE(numout,*) 'Tidal diagnostics_fast' |
|---|
| 403 | WRITE(numout,*) ' Fast Harmonic analysis ?: ln_diaharm_fast= ', ln_diaharm_fast |
|---|
| 404 | WRITE(numout,*) ' Store output in restart?: ln_diaharm_store= ', ln_diaharm_store |
|---|
| 405 | WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute |
|---|
| 406 | WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart |
|---|
| 407 | WRITE(numout,*) ' SSH harmonic analysis: ln_ana_ssh = ', ln_ana_ssh |
|---|
| 408 | WRITE(numout,*) ' Barotropic Velocities harmonic analysis: ln_ana_uvbar = ', ln_ana_uvbar |
|---|
| 409 | WRITE(numout,*) ' Bed Friction for harmonic analysis (not implemented): ln_ana_bfric = ', ln_ana_bfric |
|---|
| 410 | WRITE(numout,*) ' Density harmonic analysis: ln_ana_rho = ', ln_ana_rho |
|---|
| 411 | WRITE(numout,*) ' 3D velocities harmonic analysis: ln_ana_uv3d = ', ln_ana_uv3d |
|---|
| 412 | WRITE(numout,*) ' Vertical Velocities harmonic analysis: ln_ana_w3d = ', ln_ana_w3d |
|---|
| 413 | WRITE(numout,*) ' Names of harmonics: tname = ', tname |
|---|
| 414 | WRITE(numout,*) ' Max number of harmonics: jpmax_harmo = ', jpmax_harmo |
|---|
| 415 | WRITE(numout,*) ' Number of Harmonics: nb_harmo = ', nb_harmo |
|---|
| 416 | WRITE(numout,*) ' Multi-year harmonic analysis: ln_diaharm_multiyear = ', ln_diaharm_multiyear |
|---|
| 417 | WRITE(numout,*) ' Multi-year harmonic analysis - number of years: nn_diaharm_multiyear = ', nn_diaharm_multiyear |
|---|
| 418 | WRITE(numout,*) ' Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily |
|---|
| 419 | WRITE(numout,*) ' Number of Harmonics: nyear, nmonth = ', nyear, nmonth |
|---|
| 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 |
|---|
| 422 | |
|---|
| 423 | ENDIF |
|---|
| 424 | ! JT |
|---|
| 425 | |
|---|
| 426 | |
|---|
| 427 | IF ( ln_diaharm_multiyear ) THEN |
|---|
| 428 | ln_diaharm_store = .True. |
|---|
| 429 | IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth |
|---|
| 430 | IF ((mod(nyear,nn_diaharm_multiyear) == 0) .AND. ( nmonth == 1 )) THEN ! Jan, year = 1980,2000,2020,2040, restart tidal calculation |
|---|
| 431 | ln_diaharm_read_restart = .FALSE. |
|---|
| 432 | IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart |
|---|
| 433 | ELSE |
|---|
| 434 | ln_diaharm_read_restart = .TRUE. |
|---|
| 435 | IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart |
|---|
| 436 | ENDIF |
|---|
| 437 | |
|---|
| 438 | |
|---|
| 439 | |
|---|
| 440 | IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth |
|---|
| 441 | IF ((mod(nyear,nn_diaharm_multiyear) == (nn_diaharm_multiyear - 1)) .AND. ( nmonth == 12 )) THEN ! Dec year = 1999,2019,2039,2040, restart tidal calculation |
|---|
| 442 | ln_diaharm_compute = .TRUE. |
|---|
| 443 | IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute |
|---|
| 444 | ELSE |
|---|
| 445 | ln_diaharm_compute = .FALSE. |
|---|
| 446 | IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute |
|---|
| 447 | ENDIF |
|---|
| 448 | |
|---|
| 449 | ENDIF |
|---|
| 450 | IF ( kt < 10 ) THEN |
|---|
| 451 | ln_diaharm_read_restart = .FALSE. |
|---|
| 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 |
|---|
| 455 | ENDIF |
|---|
| 456 | |
|---|
| 457 | ! JT |
|---|
| 458 | |
|---|
| 459 | |
|---|
| 460 | ! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90 |
|---|
| 461 | nb_ana = 0 |
|---|
| 462 | DO jk=1,jpmax_harmo |
|---|
| 463 | DO ji=1,nb_harmo |
|---|
| 464 | IF(TRIM(tname(jk)) == TRIM(Wave( ntide(ji) )%cname_tide) ) THEN |
|---|
| 465 | nb_ana=nb_ana+1 |
|---|
| 466 | ENDIF |
|---|
| 467 | END DO |
|---|
| 468 | END DO |
|---|
| 469 | ! |
|---|
| 470 | IF(lwp) THEN |
|---|
| 471 | WRITE(numout,*) ' Namelist nam_diaharm_fast' |
|---|
| 472 | WRITE(numout,*) ' nb_ana = ', nb_ana |
|---|
| 473 | CALL flush(numout) |
|---|
| 474 | ENDIF |
|---|
| 475 | ! |
|---|
| 476 | IF (nb_ana > nharm_max) THEN |
|---|
| 477 | IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop' |
|---|
| 478 | IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max |
|---|
| 479 | nstop = nstop + 1 |
|---|
| 480 | ENDIF |
|---|
| 481 | |
|---|
| 482 | |
|---|
| 483 | ALLOCATE(ntide_all(nb_ana)) |
|---|
| 484 | ALLOCATE(ntide_sub(nb_ana)) |
|---|
| 485 | |
|---|
| 486 | DO jk=1,nb_ana |
|---|
| 487 | DO ji=1,nb_harmo |
|---|
| 488 | IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN |
|---|
| 489 | ntide_sub(jk) = ji |
|---|
| 490 | ntide_all(jk) = ntide(ji) |
|---|
| 491 | EXIT |
|---|
| 492 | END IF |
|---|
| 493 | END DO |
|---|
| 494 | END DO |
|---|
| 495 | |
|---|
| 496 | ALLOCATE( anau(nb_ana) ) |
|---|
| 497 | ALLOCATE( anav(nb_ana) ) |
|---|
| 498 | ALLOCATE( anaf(nb_ana) ) |
|---|
| 499 | |
|---|
| 500 | |
|---|
| 501 | IF( ln_diaharm_fast ) THEN |
|---|
| 502 | |
|---|
| 503 | ! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS |
|---|
| 504 | nvar_2d = 0; nvar_3d = 0 |
|---|
| 505 | IF ( ln_ana_ssh ) nvar_2d = nvar_2d + 1 ! analysis elevation |
|---|
| 506 | IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2 ! analysis depth-averaged velocity |
|---|
| 507 | IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2 ! analysis bottom friction |
|---|
| 508 | |
|---|
| 509 | IF ( ln_ana_rho ) nvar_3d = nvar_3d + 1 ! analysis density |
|---|
| 510 | IF ( ln_ana_uv3d ) nvar_3d = nvar_3d + 2 ! analysis 3D horizontal velocities |
|---|
| 511 | IF ( ln_ana_w3d ) nvar_3d = nvar_3d + 1 ! analysis 3D vertical velocity |
|---|
| 512 | |
|---|
| 513 | ! CHECK IF SOMETHING TO RUN |
|---|
| 514 | IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE. ! no 2d to run |
|---|
| 515 | IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE. ! no 3d to run |
|---|
| 516 | ! IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE. |
|---|
| 517 | ! IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE. |
|---|
| 518 | |
|---|
| 519 | IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN |
|---|
| 520 | |
|---|
| 521 | ! DO ALLOCATIONS |
|---|
| 522 | IF ( lk_diaharm_2D ) THEN |
|---|
| 523 | ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj, nvar_2d) ) |
|---|
| 524 | ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) ) |
|---|
| 525 | ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) ) |
|---|
| 526 | ALLOCATE( g_out2D (jpi,jpj) ) |
|---|
| 527 | ALLOCATE( h_out2D (jpi,jpj) ) |
|---|
| 528 | ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0 |
|---|
| 529 | ENDIF |
|---|
| 530 | |
|---|
| 531 | IF ( lk_diaharm_3D ) THEN |
|---|
| 532 | ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) ) |
|---|
| 533 | ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) ) |
|---|
| 534 | ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) ) |
|---|
| 535 | ALLOCATE( g_out3D (jpi,jpj,jpk) ) |
|---|
| 536 | ALLOCATE( h_out3D (jpi,jpj,jpk) ) |
|---|
| 537 | ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0 |
|---|
| 538 | ENDIF |
|---|
| 539 | |
|---|
| 540 | ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) ) |
|---|
| 541 | ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) ) |
|---|
| 542 | ALLOCATE( bzz(nb_ana*2+1) ) |
|---|
| 543 | ALLOCATE( x (nb_ana*2+1) ) |
|---|
| 544 | ALLOCATE( c (nb_ana*2+1) ) |
|---|
| 545 | ! END ALLOCATE |
|---|
| 546 | |
|---|
| 547 | ! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL |
|---|
| 548 | ! MAKES THINGS EASIER AND FASTER LATER |
|---|
| 549 | ! !!! UGLY !!! |
|---|
| 550 | jh = 1; k2d = 0; |
|---|
| 551 | IF ( ln_ana_ssh ) THEN |
|---|
| 552 | k2d = k2d + 1; m_posi_2d(k2d) = jh |
|---|
| 553 | IF(lwp) WRITE(numout,*) " - ssh harmonic analysis activated (ln_ana_ssh)" |
|---|
| 554 | ENDIF |
|---|
| 555 | jh = jh + 1 |
|---|
| 556 | IF ( ln_ana_uvbar ) THEN |
|---|
| 557 | k2d = k2d + 1; m_posi_2d(k2d) = jh |
|---|
| 558 | jh = jh + 1 |
|---|
| 559 | k2d = k2d + 1; m_posi_2d(k2d) = jh |
|---|
| 560 | IF(lwp) WRITE(numout,*) " - barotropic currents harmonic analysis activated (ln_ana_uvbar)" |
|---|
| 561 | ELSE |
|---|
| 562 | jh = jh + 1 |
|---|
| 563 | ENDIF |
|---|
| 564 | jh = jh + 1 |
|---|
| 565 | IF ( ln_ana_bfric ) THEN |
|---|
| 566 | k2d = k2d + 1; m_posi_2d(k2d) = jh |
|---|
| 567 | jh = jh + 1; |
|---|
| 568 | k2d = k2d + 1; m_posi_2d(k2d) = jh |
|---|
| 569 | IF(lwp) WRITE(numout,*) " - bottom friction harmonic analysis activated (ln_ana_vbfr)" |
|---|
| 570 | ELSE |
|---|
| 571 | jh = jh + 1 |
|---|
| 572 | ENDIF |
|---|
| 573 | |
|---|
| 574 | ! and for 3D |
|---|
| 575 | jh = 1; k3d = 0; |
|---|
| 576 | IF ( ln_ana_rho ) THEN |
|---|
| 577 | k3d = k3d + 1; m_posi_3d(k3d) = jh |
|---|
| 578 | IF(lwp) WRITE(numout,*) " - 3D density harmonic analysis activated (ln_ana_rho)" |
|---|
| 579 | ENDIF |
|---|
| 580 | jh = jh + 1 |
|---|
| 581 | IF ( ln_ana_uv3d ) THEN |
|---|
| 582 | k3d = k3d + 1; m_posi_3d(k3d) = jh |
|---|
| 583 | jh = jh + 1 |
|---|
| 584 | k3d = k3d + 1; m_posi_3d(k3d) = jh |
|---|
| 585 | IF(lwp) WRITE(numout,*) " - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)" |
|---|
| 586 | ELSE |
|---|
| 587 | jh = jh + 1 |
|---|
| 588 | ENDIF |
|---|
| 589 | jh = jh + 1 |
|---|
| 590 | IF ( ln_ana_w3d ) THEN |
|---|
| 591 | k3d = k3d + 1; m_posi_3d(k3d) = jh |
|---|
| 592 | IF(lwp) WRITE(numout,*) " - 3D vertical currents harmonic analysis activated (ln_ana_w3d)" |
|---|
| 593 | ENDIF |
|---|
| 594 | |
|---|
| 595 | ! SELECT AND STORE FREQUENCIES |
|---|
| 596 | IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency ' |
|---|
| 597 | DO jh=1,nb_ana |
|---|
| 598 | om_tide(jh) = omega_tide( ntide_sub(jh) ) |
|---|
| 599 | IF(lwp) WRITE(numout,*) ' - ',tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr" |
|---|
| 600 | ENDDO |
|---|
| 601 | |
|---|
| 602 | ! READ RESTART IF |
|---|
| 603 | IF ( ln_diaharm_read_restart ) THEN |
|---|
| 604 | IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run. kt = ",kt |
|---|
| 605 | ! Need to read in bssh bz, cc anau anav and anaf |
|---|
| 606 | call harm_rst_read ! This reads in from the previous day |
|---|
| 607 | ! Currrently the data in in assci format |
|---|
| 608 | ELSE |
|---|
| 609 | |
|---|
| 610 | IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh. kt = ",kt |
|---|
| 611 | |
|---|
| 612 | IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,: ) = 0.0_wp |
|---|
| 613 | IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp |
|---|
| 614 | cc = 0.0_wp |
|---|
| 615 | a (:,:) = 0.0_wp ! NB |
|---|
| 616 | bzz (:) = 0.0_wp |
|---|
| 617 | x (:) = 0.0_wp |
|---|
| 618 | c (:) = 0.0_wp |
|---|
| 619 | anau (:) = 0.0_wp |
|---|
| 620 | anav (:) = 0.0_wp |
|---|
| 621 | anaf (:) = 0.0_wp |
|---|
| 622 | |
|---|
| 623 | DO jh = 1, nb_ana |
|---|
| 624 | anau(jh) = utide ( ntide_sub(jh) ) |
|---|
| 625 | anav(jh) = v0tide( ntide_sub(jh) ) |
|---|
| 626 | anaf(jh) = ftide ( ntide_sub(jh) ) |
|---|
| 627 | END DO |
|---|
| 628 | |
|---|
| 629 | fjulday_startharm=fjulday !Set this at very start and store |
|---|
| 630 | !JT this is a mistake - only works on daily cycles, should use fjulnsec_dayday |
|---|
| 631 | |
|---|
| 632 | IF (lwp) THEN |
|---|
| 633 | WRITE(numout,*) '--------------------------' |
|---|
| 634 | WRITE(numout,*) ' - Output anaf for check' |
|---|
| 635 | WRITE(numout,*) 'ANA F', anaf |
|---|
| 636 | WRITE(numout,*) 'ANA U', anau |
|---|
| 637 | WRITE(numout,*) 'ANA V', anav |
|---|
| 638 | WRITE(numout,*) 'fjulday',fjulday |
|---|
| 639 | WRITE(numout,*) 'fjulday_startharm',fjulday_startharm |
|---|
| 640 | WRITE(numout,*) 'nsec_day',nsec_day |
|---|
| 641 | WRITE(numout,*) 'kt',kt |
|---|
| 642 | WRITE(numout,*) '--------------------------' |
|---|
| 643 | ENDIF |
|---|
| 644 | |
|---|
| 645 | ENDIF |
|---|
| 646 | |
|---|
| 647 | ELSE |
|---|
| 648 | |
|---|
| 649 | IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis" |
|---|
| 650 | |
|---|
| 651 | ENDIF |
|---|
| 652 | ENDIF |
|---|
| 653 | |
|---|
| 654 | END SUBROUTINE harm_ana_init |
|---|
| 655 | ! |
|---|
| 656 | SUBROUTINE gelim (a,b,x,n) |
|---|
| 657 | !!---------------------------------------------------------------------- |
|---|
| 658 | !! *** ROUTINE harm_ana *** |
|---|
| 659 | !! |
|---|
| 660 | !! ** Purpose : Guassian elimination |
|---|
| 661 | !! |
|---|
| 662 | !! |
|---|
| 663 | !! ** Action : - first action (share memory array/varible modified |
|---|
| 664 | !! in this routine |
|---|
| 665 | !! - second action ..... |
|---|
| 666 | !! - ..... |
|---|
| 667 | !! |
|---|
| 668 | !! References : |
|---|
| 669 | !! Give references if exist otherwise suppress these lines |
|---|
| 670 | !! |
|---|
| 671 | !! History : |
|---|
| 672 | implicit none |
|---|
| 673 | ! |
|---|
| 674 | integer :: n |
|---|
| 675 | REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1) |
|---|
| 676 | REAL(WP) :: x(nb_ana*2+1) |
|---|
| 677 | INTEGER :: row,col,prow,pivrow,rrow |
|---|
| 678 | REAL(WP) :: atemp |
|---|
| 679 | REAL(WP) :: pivot |
|---|
| 680 | REAL(WP) :: m |
|---|
| 681 | |
|---|
| 682 | do row=1,n-1 |
|---|
| 683 | pivrow=row |
|---|
| 684 | pivot=a(row,n-row+1) |
|---|
| 685 | do prow=row+1,n |
|---|
| 686 | if (abs(a(prow,n-row+1)).gt.abs(pivot) ) then |
|---|
| 687 | pivot=a(prow,n-row+1) |
|---|
| 688 | pivrow=prow |
|---|
| 689 | endif |
|---|
| 690 | enddo |
|---|
| 691 | ! swap row and prow |
|---|
| 692 | if ( pivrow .ne. row ) then |
|---|
| 693 | atemp=b(pivrow) |
|---|
| 694 | b(pivrow)=b(row) |
|---|
| 695 | b(row)=atemp |
|---|
| 696 | do col=1,n |
|---|
| 697 | atemp=a(pivrow,col) |
|---|
| 698 | a(pivrow,col)=a(row,col) |
|---|
| 699 | a(row,col)=atemp |
|---|
| 700 | enddo |
|---|
| 701 | endif |
|---|
| 702 | |
|---|
| 703 | do rrow=row+1,n |
|---|
| 704 | if (a(row,row).ne.0) then |
|---|
| 705 | |
|---|
| 706 | m=-a(rrow,n-row+1)/a(row,n-row+1) |
|---|
| 707 | do col=1,n |
|---|
| 708 | a(rrow,col)=m*a(row,col)+a(rrow,col) |
|---|
| 709 | enddo |
|---|
| 710 | b(rrow)=m*b(row)+b(rrow) |
|---|
| 711 | endif |
|---|
| 712 | enddo |
|---|
| 713 | enddo |
|---|
| 714 | ! back substitution now |
|---|
| 715 | |
|---|
| 716 | x(1)=b(n)/a(n,1) |
|---|
| 717 | do row=n-1,1,-1 |
|---|
| 718 | x(n-row+1)=b(row) |
|---|
| 719 | do col=1,(n-row) |
|---|
| 720 | x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) |
|---|
| 721 | enddo |
|---|
| 722 | |
|---|
| 723 | x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1)) |
|---|
| 724 | enddo |
|---|
| 725 | |
|---|
| 726 | return |
|---|
| 727 | END SUBROUTINE gelim |
|---|
| 728 | |
|---|
| 729 | SUBROUTINE harm_ana_out |
|---|
| 730 | !!---------------------------------------------------------------------- |
|---|
| 731 | !! *** ROUTINE harm_ana_init *** |
|---|
| 732 | !! |
|---|
| 733 | !! ** Purpose : initialization of .... |
|---|
| 734 | !! |
|---|
| 735 | !! ** Method : blah blah blah ... |
|---|
| 736 | !! |
|---|
| 737 | !! ** input : Namlist namexa |
|---|
| 738 | !! |
|---|
| 739 | !! ** Action : ... |
|---|
| 740 | !! |
|---|
| 741 | !! history : |
|---|
| 742 | !! 9.0 ! 03-08 (Autor Names) Original code |
|---|
| 743 | !!---------------------------------------------------------------------- |
|---|
| 744 | USE dianam ! build name of file (routine) |
|---|
| 745 | |
|---|
| 746 | !! * local declarations |
|---|
| 747 | INTEGER :: ji, jj, jk, jgrid, jh ! dummy loop indices |
|---|
| 748 | ! INTEGER :: nh_T |
|---|
| 749 | ! INTEGER :: nid_harm |
|---|
| 750 | ! CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names |
|---|
| 751 | ! CHARACTER (len=40) :: clhstnamu, clhstnamv ! temporary names |
|---|
| 752 | CHARACTER (len=40) :: suffix |
|---|
| 753 | CHARACTER (len=40) :: tmp_name |
|---|
| 754 | ! REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi ! temporary scalars |
|---|
| 755 | |
|---|
| 756 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: amp_u2d,phi_u2d, amp_v2d,phi_v2d ! arrays for output |
|---|
| 757 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:,:) :: amp_u3d,phi_u3d, amp_v3d,phi_v3d ! arrays for output |
|---|
| 758 | |
|---|
| 759 | REAL(wp) :: tmp_u_amp ,tmp_v_amp ,tmp_u_phi ,tmp_v_phi |
|---|
| 760 | REAL(wp) :: a_u, b_u, a_v, b_v, twodelta, delta, alpha2, alpha, qmin, qmax, ecc,thetamax, thetamin |
|---|
| 761 | REAL(wp) :: Qc, Qac, gc,gac, Phi_Ua, dir_Ua, polarity |
|---|
| 762 | REAL(wp) :: u_off, v_off |
|---|
| 763 | REAL(wp) :: tmpreal |
|---|
| 764 | |
|---|
| 765 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: tmp_u_amp_2d_mat,tmp_v_amp_2d_mat,tmp_u_phi_2d_mat,tmp_v_phi_2d_mat |
|---|
| 766 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: TA_u_off_t_uvbar, TA_v_off_t_uvbar |
|---|
| 767 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: TA_u_off_uvbar, TA_v_off_uvbar |
|---|
| 768 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: a_u_2d_mat,b_u_2d_mat,a_v_2d_mat,b_v_2d_mat |
|---|
| 769 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: qmax_2d_mat,qmin_2d_mat,ecc_2d_mat |
|---|
| 770 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: thetamax_2d_mat,thetamin_2d_mat,Qc_2d_mat,Qac_2d_mat |
|---|
| 771 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: gc_2d_mat,gac_2d_mat,Phi_Ua_2d_mat,dir_Ua_2d_mat |
|---|
| 772 | REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: polarity_2d_mat |
|---|
| 773 | |
|---|
| 774 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: tmp_u_amp_3d_mat,tmp_v_amp_3d_mat,tmp_u_phi_3d_mat,tmp_v_phi_3d_mat |
|---|
| 775 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: TA_u_off_t_uv3d, TA_v_off_t_uv3d |
|---|
| 776 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: TA_u_off_uv3d, TA_v_off_uv3d |
|---|
| 777 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: a_u_3d_mat,b_u_3d_mat,a_v_3d_mat,b_v_3d_mat |
|---|
| 778 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: qmax_3d_mat,qmin_3d_mat,ecc_3d_mat |
|---|
| 779 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: thetamax_3d_mat,thetamin_3d_mat,Qc_3d_mat,Qac_3d_mat |
|---|
| 780 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: gc_3d_mat,gac_3d_mat,Phi_Ua_3d_mat,dir_Ua_3d_mat |
|---|
| 781 | REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: polarity_3d_mat |
|---|
| 782 | |
|---|
| 783 | |
|---|
| 784 | |
|---|
| 785 | |
|---|
| 786 | IF (ln_diaharm_postproc_vel) THEN |
|---|
| 787 | IF (ln_ana_uvbar) THEN |
|---|
| 788 | 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) ) |
|---|
| 789 | |
|---|
| 790 | ALLOCATE( TA_u_off_t_uvbar(jpi,jpj), TA_v_off_t_uvbar(jpi,jpj) ) |
|---|
| 791 | ALLOCATE( TA_u_off_uvbar(jpi,jpj), TA_v_off_uvbar(jpi,jpj) ) |
|---|
| 792 | |
|---|
| 793 | 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)) |
|---|
| 794 | 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)) |
|---|
| 795 | ALLOCATE(qmax_2d_mat(jpi,jpj),qmin_2d_mat(jpi,jpj),ecc_2d_mat(jpi,jpj)) |
|---|
| 796 | ALLOCATE(thetamax_2d_mat(jpi,jpj),thetamin_2d_mat(jpi,jpj),Qc_2d_mat(jpi,jpj),Qac_2d_mat(jpi,jpj)) |
|---|
| 797 | ALLOCATE(gc_2d_mat(jpi,jpj),gac_2d_mat(jpi,jpj),Phi_Ua_2d_mat(jpi,jpj),dir_Ua_2d_mat(jpi,jpj)) |
|---|
| 798 | ALLOCATE(polarity_2d_mat(jpi,jpj)) |
|---|
| 799 | |
|---|
| 800 | ENDIF |
|---|
| 801 | |
|---|
| 802 | |
|---|
| 803 | IF (ln_ana_uv3d) THEN |
|---|
| 804 | 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) ) |
|---|
| 805 | |
|---|
| 806 | ALLOCATE( TA_u_off_t_uv3d(jpi,jpj,jpk), TA_v_off_t_uv3d(jpi,jpj,jpk) ) |
|---|
| 807 | ALLOCATE( TA_u_off_uv3d(jpi,jpj,jpk), TA_v_off_uv3d(jpi,jpj,jpk) ) |
|---|
| 808 | |
|---|
| 809 | 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)) |
|---|
| 810 | 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)) |
|---|
| 811 | ALLOCATE(qmax_3d_mat(jpi,jpj,jpk),qmin_3d_mat(jpi,jpj,jpk),ecc_3d_mat(jpi,jpj,jpk)) |
|---|
| 812 | ALLOCATE(thetamax_3d_mat(jpi,jpj,jpk),thetamin_3d_mat(jpi,jpj,jpk),Qc_3d_mat(jpi,jpj,jpk),Qac_3d_mat(jpi,jpj,jpk)) |
|---|
| 813 | ALLOCATE(gc_3d_mat(jpi,jpj,jpk),gac_3d_mat(jpi,jpj,jpk),Phi_Ua_3d_mat(jpi,jpj,jpk),dir_Ua_3d_mat(jpi,jpj,jpk)) |
|---|
| 814 | ALLOCATE(polarity_3d_mat(jpi,jpj,jpk)) |
|---|
| 815 | |
|---|
| 816 | ENDIF |
|---|
| 817 | ENDIF |
|---|
| 818 | |
|---|
| 819 | do jgrid=1,nvar_2d |
|---|
| 820 | do jh=1,nb_ana |
|---|
| 821 | h_out2D = 0.0 |
|---|
| 822 | g_out2D = 0.0 |
|---|
| 823 | do jj=1,nlcj |
|---|
| 824 | do ji=1,nlci |
|---|
| 825 | cca=g_cosamp2D(jh,ji,jj,jgrid) |
|---|
| 826 | ssa=g_sinamp2D(jh,ji,jj,jgrid) |
|---|
| 827 | h_out2D(ji,jj)=sqrt(cca**2+ssa**2) |
|---|
| 828 | IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN |
|---|
| 829 | g_out2D(ji,jj)= 0.0_wp |
|---|
| 830 | ELSE |
|---|
| 831 | g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca) |
|---|
| 832 | ENDIF |
|---|
| 833 | IF (h_out2D(ji,jj).ne.0) THEN |
|---|
| 834 | h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh) |
|---|
| 835 | ENDIF |
|---|
| 836 | IF (g_out2D(ji,jj).ne.0) THEN !Correct and take modulus |
|---|
| 837 | !JT |
|---|
| 838 | !JT g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0) |
|---|
| 839 | !JT |
|---|
| 840 | g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh))/rad , 360.0) |
|---|
| 841 | if (g_out2D(ji,jj).gt.360.0) then |
|---|
| 842 | g_out2D(ji,jj)=g_out2D(ji,jj)-360.0 |
|---|
| 843 | else if (g_out2D(ji,jj).lt.0.0) then |
|---|
| 844 | g_out2D(ji,jj)=g_out2D(ji,jj)+360.0 |
|---|
| 845 | endif |
|---|
| 846 | ENDIF |
|---|
| 847 | enddo |
|---|
| 848 | enddo |
|---|
| 849 | ! |
|---|
| 850 | ! NETCDF OUTPUT |
|---|
| 851 | suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) ) |
|---|
| 852 | IF(lwp) WRITE(numout,*) "diaharm_fast", suffix |
|---|
| 853 | |
|---|
| 854 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) |
|---|
| 855 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 856 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D) |
|---|
| 857 | 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" |
|---|
| 858 | CALL iom_put( TRIM(tmp_name), h_out2D(:,:) ) |
|---|
| 859 | ELSE |
|---|
| 860 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 861 | ENDIF |
|---|
| 862 | |
|---|
| 863 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) |
|---|
| 864 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 865 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D) |
|---|
| 866 | CALL iom_put( TRIM(tmp_name), g_out2D(:,:) ) |
|---|
| 867 | ELSE |
|---|
| 868 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 869 | ENDIF |
|---|
| 870 | |
|---|
| 871 | |
|---|
| 872 | |
|---|
| 873 | IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN |
|---|
| 874 | |
|---|
| 875 | !IF (m_posi_2d(jgrid) == 2) THEN |
|---|
| 876 | IF (TRIM(suffix) == TRIM('u2d')) THEN |
|---|
| 877 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u2d '//TRIM(suffix) |
|---|
| 878 | do jj=1,nlcj |
|---|
| 879 | do ji=1,nlci |
|---|
| 880 | if (ssumask(ji,jj) == 1) THEN |
|---|
| 881 | amp_u2d(jh,ji,jj) = h_out2D(ji,jj) |
|---|
| 882 | phi_u2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0 |
|---|
| 883 | else |
|---|
| 884 | amp_u2d(jh,ji,jj) = 0. |
|---|
| 885 | phi_u2d(jh,ji,jj) = 0. |
|---|
| 886 | ENDIF |
|---|
| 887 | enddo |
|---|
| 888 | enddo |
|---|
| 889 | ENDIF |
|---|
| 890 | |
|---|
| 891 | !IF (m_posi_2d(jgrid) == 3) THEN |
|---|
| 892 | IF (TRIM(suffix) == TRIM('v2d')) THEN |
|---|
| 893 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v2d '//TRIM(suffix) |
|---|
| 894 | do jj=1,nlcj |
|---|
| 895 | do ji=1,nlci |
|---|
| 896 | if (ssvmask(ji,jj) == 1) THEN |
|---|
| 897 | amp_v2d(jh,ji,jj) = h_out2D(ji,jj) |
|---|
| 898 | phi_v2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0 |
|---|
| 899 | else |
|---|
| 900 | amp_v2d(jh,ji,jj) = 0 |
|---|
| 901 | phi_v2d(jh,ji,jj) = 0 |
|---|
| 902 | ENDIF |
|---|
| 903 | enddo |
|---|
| 904 | enddo |
|---|
| 905 | ENDIF |
|---|
| 906 | ENDIF |
|---|
| 907 | |
|---|
| 908 | CALL FLUSH(numout) |
|---|
| 909 | |
|---|
| 910 | |
|---|
| 911 | enddo |
|---|
| 912 | |
|---|
| 913 | suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) ) |
|---|
| 914 | tmp_name='TA_'//TRIM(suffix)//'_off' |
|---|
| 915 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 916 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 917 | CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid)) |
|---|
| 918 | ELSE |
|---|
| 919 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 920 | ENDIF |
|---|
| 921 | |
|---|
| 922 | CALL FLUSH(numout) |
|---|
| 923 | |
|---|
| 924 | |
|---|
| 925 | |
|---|
| 926 | |
|---|
| 927 | |
|---|
| 928 | |
|---|
| 929 | IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN |
|---|
| 930 | |
|---|
| 931 | !IF (m_posi_2d(jgrid) == 2) THEN |
|---|
| 932 | IF (TRIM(suffix) == TRIM('u2d')) THEN |
|---|
| 933 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: TA_u_off_uvbar" |
|---|
| 934 | do jj=1,nlcj |
|---|
| 935 | do ji=1,nlci |
|---|
| 936 | if (ssumask(ji,jj) == 1) THEN |
|---|
| 937 | TA_u_off_uvbar(ji,jj) = g_cosamp2D( 0,ji,jj,jgrid) |
|---|
| 938 | else |
|---|
| 939 | TA_u_off_uvbar(ji,jj) = 0. |
|---|
| 940 | ENDIF |
|---|
| 941 | enddo !ji |
|---|
| 942 | enddo !jj |
|---|
| 943 | ENDIF !u2d |
|---|
| 944 | |
|---|
| 945 | !IF (m_posi_2d(jgrid) == 3) THEN |
|---|
| 946 | IF (TRIM(suffix) == TRIM('v2d')) THEN |
|---|
| 947 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: TA_v_off_uvbar" |
|---|
| 948 | do jj=1,nlcj |
|---|
| 949 | do ji=1,nlci |
|---|
| 950 | if (ssvmask(ji,jj) == 1) THEN |
|---|
| 951 | TA_v_off_uvbar(ji,jj) = g_cosamp2D( 0,ji,jj,jgrid) |
|---|
| 952 | else |
|---|
| 953 | TA_v_off_uvbar(ji,jj) = 0. |
|---|
| 954 | ENDIF |
|---|
| 955 | enddo !ji |
|---|
| 956 | enddo !jj |
|---|
| 957 | ENDIF !uvd |
|---|
| 958 | |
|---|
| 959 | |
|---|
| 960 | ENDIF ! ln_diaharm_postproc_vel .AND. ln_ana_uvbar |
|---|
| 961 | |
|---|
| 962 | |
|---|
| 963 | enddo ! jgrid=1,nvar_2d |
|---|
| 964 | ! |
|---|
| 965 | ! DO THE SAME FOR 3D VARIABLES |
|---|
| 966 | ! |
|---|
| 967 | do jgrid=1,nvar_3d |
|---|
| 968 | do jh=1,nb_ana |
|---|
| 969 | h_out3D = 0.0 |
|---|
| 970 | g_out3D = 0.0 |
|---|
| 971 | DO jk=1,jpkm1 |
|---|
| 972 | do jj=1,nlcj |
|---|
| 973 | do ji=1,nlci |
|---|
| 974 | cca=g_cosamp3D(jh,ji,jj,jk,jgrid) |
|---|
| 975 | ssa=g_sinamp3D(jh,ji,jj,jk,jgrid) |
|---|
| 976 | h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2) |
|---|
| 977 | IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN |
|---|
| 978 | g_out3D(ji,jj,jk) = 0.0_wp |
|---|
| 979 | ELSE |
|---|
| 980 | g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca) |
|---|
| 981 | ENDIF |
|---|
| 982 | IF (h_out3D(ji,jj,jk).ne.0) THEN |
|---|
| 983 | h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh) |
|---|
| 984 | ENDIF |
|---|
| 985 | IF (g_out3D(ji,jj,jk).ne.0) THEN !Correct and take modulus |
|---|
| 986 | !JT |
|---|
| 987 | !JT g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0) |
|---|
| 988 | !JT |
|---|
| 989 | g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh))/rad , 360.0) |
|---|
| 990 | if (g_out3D(ji,jj,jk).gt.360.0) then |
|---|
| 991 | g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0 |
|---|
| 992 | else if (g_out3D(ji,jj,jk).lt.0.0) then |
|---|
| 993 | g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0 |
|---|
| 994 | endif |
|---|
| 995 | ENDIF |
|---|
| 996 | enddo ! ji |
|---|
| 997 | enddo ! jj |
|---|
| 998 | ENDDO ! jk |
|---|
| 999 | ! |
|---|
| 1000 | ! NETCDF OUTPUT |
|---|
| 1001 | suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) ) |
|---|
| 1002 | IF(lwp) WRITE(numout,*) "diaharm_fast", suffix |
|---|
| 1003 | |
|---|
| 1004 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix) |
|---|
| 1005 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1006 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D) |
|---|
| 1007 | CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) ) |
|---|
| 1008 | ELSE |
|---|
| 1009 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 1010 | ENDIF |
|---|
| 1011 | |
|---|
| 1012 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix) |
|---|
| 1013 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1014 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D) |
|---|
| 1015 | CALL iom_put(tmp_name, g_out3D(:,:,:) ) |
|---|
| 1016 | ELSE |
|---|
| 1017 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 1018 | ENDIF |
|---|
| 1019 | |
|---|
| 1020 | |
|---|
| 1021 | |
|---|
| 1022 | IF (ln_diaharm_postproc_vel .AND. ln_ana_uv3d) THEN |
|---|
| 1023 | |
|---|
| 1024 | !IF (m_posi_2d(jgrid) == 2) THEN |
|---|
| 1025 | IF (TRIM(suffix) == TRIM('u3d')) THEN |
|---|
| 1026 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u3d '//TRIM(suffix) |
|---|
| 1027 | DO jk=1,jpkm1 |
|---|
| 1028 | do jj=1,nlcj |
|---|
| 1029 | do ji=1,nlci |
|---|
| 1030 | if (umask(ji,jj,jk) == 1) THEN |
|---|
| 1031 | amp_u3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) |
|---|
| 1032 | phi_u3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0 |
|---|
| 1033 | else |
|---|
| 1034 | amp_u3d(jh,ji,jj,jk) = 0. |
|---|
| 1035 | phi_u3d(jh,ji,jj,jk) = 0. |
|---|
| 1036 | ENDIF |
|---|
| 1037 | enddo |
|---|
| 1038 | enddo |
|---|
| 1039 | enddo |
|---|
| 1040 | ENDIF |
|---|
| 1041 | |
|---|
| 1042 | !IF (m_posi_2d(jgrid) == 3) THEN |
|---|
| 1043 | IF (TRIM(suffix) == TRIM('v3d')) THEN |
|---|
| 1044 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v3d '//TRIM(suffix) |
|---|
| 1045 | DO jk=1,jpkm1 |
|---|
| 1046 | do jj=1,nlcj |
|---|
| 1047 | do ji=1,nlci |
|---|
| 1048 | if (vmask(ji,jj,jk) == 1) THEN |
|---|
| 1049 | amp_v3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) |
|---|
| 1050 | phi_v3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0 |
|---|
| 1051 | else |
|---|
| 1052 | amp_v3d(jh,ji,jj,jk) = 0. |
|---|
| 1053 | phi_v3d(jh,ji,jj,jk) = 0. |
|---|
| 1054 | ENDIF |
|---|
| 1055 | enddo |
|---|
| 1056 | enddo |
|---|
| 1057 | enddo |
|---|
| 1058 | ENDIF |
|---|
| 1059 | ENDIF |
|---|
| 1060 | |
|---|
| 1061 | CALL FLUSH(numout) |
|---|
| 1062 | |
|---|
| 1063 | |
|---|
| 1064 | enddo ! jh |
|---|
| 1065 | |
|---|
| 1066 | suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) ) |
|---|
| 1067 | tmp_name='TA_'//TRIM(suffix)//'_off' |
|---|
| 1068 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1069 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1070 | CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid)) |
|---|
| 1071 | ELSE |
|---|
| 1072 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name) |
|---|
| 1073 | ENDIF |
|---|
| 1074 | |
|---|
| 1075 | |
|---|
| 1076 | |
|---|
| 1077 | |
|---|
| 1078 | |
|---|
| 1079 | |
|---|
| 1080 | |
|---|
| 1081 | IF (ln_diaharm_postproc_vel .AND. ln_ana_uv3d) THEN |
|---|
| 1082 | |
|---|
| 1083 | !IF (m_posi_2d(jgrid) == 2) THEN |
|---|
| 1084 | IF (TRIM(suffix) == TRIM('u3d')) THEN |
|---|
| 1085 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: TA_u_off_uv3d" |
|---|
| 1086 | DO jk=1,jpkm1 |
|---|
| 1087 | do jj=1,nlcj |
|---|
| 1088 | do ji=1,nlci |
|---|
| 1089 | if (umask(ji,jj,jk) == 1) THEN |
|---|
| 1090 | TA_u_off_uv3d(ji,jj,jk) = g_cosamp3D( 0,ji,jj,jk,jgrid) |
|---|
| 1091 | else |
|---|
| 1092 | TA_u_off_uv3d(ji,jj,jk) = 0. |
|---|
| 1093 | ENDIF |
|---|
| 1094 | enddo |
|---|
| 1095 | enddo |
|---|
| 1096 | enddo |
|---|
| 1097 | ENDIF |
|---|
| 1098 | |
|---|
| 1099 | !IF (m_posi_2d(jgrid) == 3) THEN |
|---|
| 1100 | IF (TRIM(suffix) == TRIM('v3d')) THEN |
|---|
| 1101 | if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: TA_v_off_uv3d" |
|---|
| 1102 | DO jk=1,jpkm1 |
|---|
| 1103 | do jj=1,nlcj |
|---|
| 1104 | do ji=1,nlci |
|---|
| 1105 | if (vmask(ji,jj,jk) == 1) THEN |
|---|
| 1106 | TA_v_off_uv3d(ji,jj,jk) = g_cosamp3D( 0,ji,jj,jk,jgrid) |
|---|
| 1107 | else |
|---|
| 1108 | TA_v_off_uv3d(ji,jj,jk) = 0. |
|---|
| 1109 | ENDIF |
|---|
| 1110 | enddo !jk |
|---|
| 1111 | enddo !ji |
|---|
| 1112 | enddo !jj |
|---|
| 1113 | ENDIF !uvd |
|---|
| 1114 | |
|---|
| 1115 | |
|---|
| 1116 | ENDIF ! ln_diaharm_postproc_vel .AND. ln_ana_uv3d |
|---|
| 1117 | |
|---|
| 1118 | |
|---|
| 1119 | enddo ! jgrid=1,nvar_2d |
|---|
| 1120 | |
|---|
| 1121 | |
|---|
| 1122 | |
|---|
| 1123 | CALL FLUSH(numout) |
|---|
| 1124 | |
|---|
| 1125 | |
|---|
| 1126 | |
|---|
| 1127 | |
|---|
| 1128 | IF (ln_diaharm_postproc_vel ) THEN |
|---|
| 1129 | |
|---|
| 1130 | |
|---|
| 1131 | TA_u_off_t_uvbar(:,:) = 0. |
|---|
| 1132 | TA_v_off_t_uvbar(:,:) = 0. |
|---|
| 1133 | |
|---|
| 1134 | DO jj = 1, nlcj !- 1 |
|---|
| 1135 | DO ji = 1, nlci ! - 1 |
|---|
| 1136 | |
|---|
| 1137 | IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN |
|---|
| 1138 | |
|---|
| 1139 | if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then |
|---|
| 1140 | u_off = ((TA_u_off_uvbar(ji,jj)*ssumask(ji,jj)) + (TA_u_off_uvbar(ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj)) |
|---|
| 1141 | else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then |
|---|
| 1142 | u_off = (TA_u_off_uvbar(ji,jj)*ssumask(ji,jj)) |
|---|
| 1143 | else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then |
|---|
| 1144 | u_off = (TA_u_off_uvbar(ji-1,jj)*ssumask(ji-1,jj)) |
|---|
| 1145 | else |
|---|
| 1146 | cycle |
|---|
| 1147 | end if |
|---|
| 1148 | |
|---|
| 1149 | |
|---|
| 1150 | if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then |
|---|
| 1151 | v_off = ((TA_v_off_uvbar(ji,jj)*ssvmask(ji,jj)) + (TA_v_off_uvbar(ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1)) |
|---|
| 1152 | else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then |
|---|
| 1153 | v_off = (TA_v_off_uvbar(ji,jj)*ssvmask(ji,jj)) |
|---|
| 1154 | else if ( (ssvmask(ji,jj) == 0) .AND. (ssvmask(ji,jj-1) == 1)) then |
|---|
| 1155 | v_off = (TA_v_off_uvbar(ji,jj-1)*ssvmask(ji,jj-1)) |
|---|
| 1156 | else |
|---|
| 1157 | cycle |
|---|
| 1158 | end if |
|---|
| 1159 | |
|---|
| 1160 | TA_u_off_t_uvbar(ji,jj) = u_off |
|---|
| 1161 | TA_v_off_t_uvbar(ji,jj) = v_off |
|---|
| 1162 | |
|---|
| 1163 | ENDIF |
|---|
| 1164 | END DO !ji |
|---|
| 1165 | END DO !jj |
|---|
| 1166 | |
|---|
| 1167 | tmp_name='TA_u_off_t_uvbar' |
|---|
| 1168 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1169 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1170 | CALL iom_put( TRIM(tmp_name), TA_u_off_t_uvbar(:,:)) |
|---|
| 1171 | ENDIF |
|---|
| 1172 | tmp_name='TA_v_off_t_uvbar' |
|---|
| 1173 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1174 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1175 | CALL iom_put( TRIM(tmp_name), TA_v_off_t_uvbar(:,:)) |
|---|
| 1176 | ENDIF |
|---|
| 1177 | |
|---|
| 1178 | TA_u_off_t_uvbar(:,:) = 0. |
|---|
| 1179 | TA_v_off_t_uvbar(:,:) = 0. |
|---|
| 1180 | |
|---|
| 1181 | |
|---|
| 1182 | |
|---|
| 1183 | TA_u_off_t_uv3d(:,:,:) = 0. |
|---|
| 1184 | TA_v_off_t_uv3d(:,:,:) = 0. |
|---|
| 1185 | |
|---|
| 1186 | DO jk=1,jpkm1 |
|---|
| 1187 | DO jj = 1, nlcj !- 1 |
|---|
| 1188 | DO ji = 1, nlci ! - 1 |
|---|
| 1189 | |
|---|
| 1190 | IF ( ((umask(ji,jj,jk) + umask(ji-1,jj,jk)) > 0 ) .AND. ((vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) > 0 ) ) THEN |
|---|
| 1191 | |
|---|
| 1192 | if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then |
|---|
| 1193 | u_off = ((TA_u_off_uv3d(ji,jj,jk)*umask(ji,jj,jk)) + (TA_u_off_uv3d(ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk)) |
|---|
| 1194 | else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then |
|---|
| 1195 | u_off = (TA_u_off_uv3d(ji,jj,jk)*umask(ji,jj,jk)) |
|---|
| 1196 | else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then |
|---|
| 1197 | u_off = (TA_u_off_uv3d(ji-1,jj,jk)*umask(ji-1,jj,jk)) |
|---|
| 1198 | else |
|---|
| 1199 | cycle |
|---|
| 1200 | end if |
|---|
| 1201 | |
|---|
| 1202 | |
|---|
| 1203 | if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then |
|---|
| 1204 | v_off = ((TA_v_off_uv3d(ji,jj,jk)*vmask(ji,jj,jk)) + (TA_v_off_uv3d(ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) |
|---|
| 1205 | else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then |
|---|
| 1206 | v_off = (TA_v_off_uv3d(ji,jj,jk)*vmask(ji,jj,jk)) |
|---|
| 1207 | else if ( (vmask(ji,jj,jk) == 0) .AND. (vmask(ji,jj-1,jk) == 1)) then |
|---|
| 1208 | v_off = (TA_v_off_uv3d(ji,jj-1,jk)*vmask(ji,jj-1,jk)) |
|---|
| 1209 | else |
|---|
| 1210 | cycle |
|---|
| 1211 | end if |
|---|
| 1212 | |
|---|
| 1213 | TA_u_off_t_uv3d(ji,jj,jk) = u_off |
|---|
| 1214 | TA_v_off_t_uv3d(ji,jj,jk) = v_off |
|---|
| 1215 | |
|---|
| 1216 | ENDIF |
|---|
| 1217 | END DO !ji |
|---|
| 1218 | END DO !jj |
|---|
| 1219 | END DO !jk |
|---|
| 1220 | |
|---|
| 1221 | tmp_name='TA_u_off_t_uv3d' |
|---|
| 1222 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1223 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1224 | CALL iom_put( TRIM(tmp_name), TA_u_off_t_uv3d(:,:,:)) |
|---|
| 1225 | ENDIF |
|---|
| 1226 | tmp_name='TA_v_off_t_uv3d' |
|---|
| 1227 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1228 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1229 | CALL iom_put( TRIM(tmp_name), TA_v_off_t_uv3d(:,:,:)) |
|---|
| 1230 | ENDIF |
|---|
| 1231 | |
|---|
| 1232 | TA_u_off_t_uv3d(:,:,:) = 0. |
|---|
| 1233 | TA_v_off_t_uv3d(:,:,:) = 0. |
|---|
| 1234 | |
|---|
| 1235 | ENDIF |
|---|
| 1236 | |
|---|
| 1237 | IF (ln_diaharm_postproc_vel ) THEN |
|---|
| 1238 | IF ( ln_ana_uvbar) THEN |
|---|
| 1239 | IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess barotropic velocity tidal parameters" |
|---|
| 1240 | CALL FLUSH(numout) |
|---|
| 1241 | DO jh=1,nb_ana |
|---|
| 1242 | |
|---|
| 1243 | |
|---|
| 1244 | tmp_u_amp_2d_mat(:,:) = 0. |
|---|
| 1245 | tmp_v_amp_2d_mat(:,:) = 0. |
|---|
| 1246 | tmp_u_phi_2d_mat(:,:) = 0. |
|---|
| 1247 | tmp_v_phi_2d_mat(:,:) = 0. |
|---|
| 1248 | |
|---|
| 1249 | a_u_2d_mat(:,:) = 0. |
|---|
| 1250 | b_u_2d_mat(:,:) = 0. |
|---|
| 1251 | a_v_2d_mat(:,:) = 0. |
|---|
| 1252 | b_v_2d_mat(:,:) = 0. |
|---|
| 1253 | |
|---|
| 1254 | qmax_2d_mat(:,:) = 0. |
|---|
| 1255 | qmin_2d_mat(:,:) = 0. |
|---|
| 1256 | |
|---|
| 1257 | ecc_2d_mat(:,:) = 0. |
|---|
| 1258 | thetamax_2d_mat(:,:) =0. |
|---|
| 1259 | thetamin_2d_mat(:,:) = 0. |
|---|
| 1260 | |
|---|
| 1261 | Qc_2d_mat(:,:) = 0. |
|---|
| 1262 | Qac_2d_mat(:,:) = 0. |
|---|
| 1263 | gc_2d_mat(:,:) = 0. |
|---|
| 1264 | gac_2d_mat(:,:) = 0. |
|---|
| 1265 | |
|---|
| 1266 | Phi_Ua_2d_mat(:,:) = 0. |
|---|
| 1267 | dir_Ua_2d_mat(:,:) = 0. |
|---|
| 1268 | polarity_2d_mat(:,:) = 0. |
|---|
| 1269 | |
|---|
| 1270 | |
|---|
| 1271 | DO jj = 1, nlcj !- 1 |
|---|
| 1272 | DO ji = 1, nlci ! - 1 |
|---|
| 1273 | |
|---|
| 1274 | IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN |
|---|
| 1275 | |
|---|
| 1276 | if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then |
|---|
| 1277 | |
|---|
| 1278 | 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)) |
|---|
| 1279 | 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)))) |
|---|
| 1280 | else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then |
|---|
| 1281 | tmp_u_amp = (amp_u2d(jh,ji,jj)*ssumask(ji,jj)) |
|---|
| 1282 | tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj)) |
|---|
| 1283 | else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then |
|---|
| 1284 | tmp_u_amp = (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) |
|---|
| 1285 | tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)) |
|---|
| 1286 | else |
|---|
| 1287 | cycle |
|---|
| 1288 | end if |
|---|
| 1289 | |
|---|
| 1290 | |
|---|
| 1291 | if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then |
|---|
| 1292 | 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)) |
|---|
| 1293 | 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)))) |
|---|
| 1294 | else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then |
|---|
| 1295 | tmp_v_amp = (amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) |
|---|
| 1296 | tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) |
|---|
| 1297 | else if ( (ssvmask(ji,jj) == 0) .AND. (ssvmask(ji,jj-1) == 1)) then |
|---|
| 1298 | tmp_v_amp = (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)) |
|---|
| 1299 | tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)) |
|---|
| 1300 | else |
|---|
| 1301 | cycle |
|---|
| 1302 | end if |
|---|
| 1303 | |
|---|
| 1304 | |
|---|
| 1305 | a_u = tmp_U_amp * cos(tmp_U_phi) |
|---|
| 1306 | b_u = tmp_U_amp * sin(tmp_U_phi) |
|---|
| 1307 | a_v = tmp_V_amp * cos(tmp_V_phi) |
|---|
| 1308 | b_v = tmp_V_amp * sin(tmp_V_phi) |
|---|
| 1309 | |
|---|
| 1310 | 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) ) ) ) |
|---|
| 1311 | delta = twodelta/2. |
|---|
| 1312 | |
|---|
| 1313 | !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)) ) |
|---|
| 1314 | 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)) |
|---|
| 1315 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1316 | |
|---|
| 1317 | !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)) ) |
|---|
| 1318 | alpha2 = sqrt( tmpreal ) |
|---|
| 1319 | if (alpha2 < 0) alpha2 = 0 !CYCLE |
|---|
| 1320 | alpha= sqrt( alpha2 ) |
|---|
| 1321 | |
|---|
| 1322 | |
|---|
| 1323 | !major and minor axis of the ellipse |
|---|
| 1324 | qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) |
|---|
| 1325 | |
|---|
| 1326 | tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 |
|---|
| 1327 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1328 | !qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. |
|---|
| 1329 | qmin = sqrt( tmpreal ) ! but always positive. |
|---|
| 1330 | |
|---|
| 1331 | !eccentricity of ellipse |
|---|
| 1332 | tmpreal = (qmax + qmin) |
|---|
| 1333 | if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE |
|---|
| 1334 | !ecc = (qmax - qmin)/(qmax + qmin) |
|---|
| 1335 | ecc = (qmax - qmin)/(tmpreal) |
|---|
| 1336 | |
|---|
| 1337 | ! Angle of major and minor ellipse |
|---|
| 1338 | thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) |
|---|
| 1339 | thetamin = thetamax + rpi/2. |
|---|
| 1340 | |
|---|
| 1341 | |
|---|
| 1342 | |
|---|
| 1343 | ! Rotary current components: Pugh A3.10 |
|---|
| 1344 | ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors |
|---|
| 1345 | ! so Qc = clockwise = anticyclonic = negative |
|---|
| 1346 | ! and Qac = anticlockwise = cyclonic = negative |
|---|
| 1347 | |
|---|
| 1348 | tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) |
|---|
| 1349 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1350 | !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)) ) |
|---|
| 1351 | Qc = 0.5*sqrt( tmpreal ) |
|---|
| 1352 | |
|---|
| 1353 | tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) |
|---|
| 1354 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1355 | !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)) ) |
|---|
| 1356 | Qac = 0.5*sqrt( tmpreal ) |
|---|
| 1357 | |
|---|
| 1358 | |
|---|
| 1359 | gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & |
|---|
| 1360 | & ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) |
|---|
| 1361 | gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & |
|---|
| 1362 | & ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) |
|---|
| 1363 | |
|---|
| 1364 | !Pugh A3.2 |
|---|
| 1365 | Phi_Ua = -0.5*(gac - gc) |
|---|
| 1366 | dir_Ua = 0.5*(gac + gc) ! positive from x axis |
|---|
| 1367 | |
|---|
| 1368 | tmpreal = qmax |
|---|
| 1369 | !if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE |
|---|
| 1370 | if (tmpreal == 0) then |
|---|
| 1371 | polarity = 0 |
|---|
| 1372 | else |
|---|
| 1373 | polarity = (Qac - Qc)/qmax |
|---|
| 1374 | endif |
|---|
| 1375 | |
|---|
| 1376 | tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp |
|---|
| 1377 | tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp |
|---|
| 1378 | tmp_u_phi_2d_mat(ji,jj) = tmp_u_phi |
|---|
| 1379 | tmp_v_phi_2d_mat(ji,jj) = tmp_v_phi |
|---|
| 1380 | |
|---|
| 1381 | |
|---|
| 1382 | a_u_2d_mat(ji,jj) = a_u |
|---|
| 1383 | b_u_2d_mat(ji,jj) = b_u |
|---|
| 1384 | a_v_2d_mat(ji,jj) = a_v |
|---|
| 1385 | b_v_2d_mat(ji,jj) = b_v |
|---|
| 1386 | |
|---|
| 1387 | qmax_2d_mat(ji,jj) = qmax |
|---|
| 1388 | qmin_2d_mat(ji,jj) = qmin |
|---|
| 1389 | |
|---|
| 1390 | ecc_2d_mat(ji,jj) = ecc |
|---|
| 1391 | thetamax_2d_mat(ji,jj) = thetamax |
|---|
| 1392 | thetamin_2d_mat(ji,jj) = thetamin |
|---|
| 1393 | |
|---|
| 1394 | Qc_2d_mat(ji,jj) = Qc |
|---|
| 1395 | Qac_2d_mat(ji,jj) = Qac |
|---|
| 1396 | gc_2d_mat(ji,jj) = gc |
|---|
| 1397 | gac_2d_mat(ji,jj) = gac |
|---|
| 1398 | |
|---|
| 1399 | Phi_Ua_2d_mat(ji,jj) = Phi_Ua |
|---|
| 1400 | dir_Ua_2d_mat(ji,jj) = dir_Ua |
|---|
| 1401 | polarity_2d_mat(ji,jj) = polarity |
|---|
| 1402 | |
|---|
| 1403 | ENDIF |
|---|
| 1404 | END DO |
|---|
| 1405 | END DO |
|---|
| 1406 | |
|---|
| 1407 | |
|---|
| 1408 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar' |
|---|
| 1409 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1410 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1411 | CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:)) |
|---|
| 1412 | ENDIF |
|---|
| 1413 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar' |
|---|
| 1414 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1415 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1416 | CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:)) |
|---|
| 1417 | ENDIF |
|---|
| 1418 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar' |
|---|
| 1419 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1420 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1421 | CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:)) |
|---|
| 1422 | ENDIF |
|---|
| 1423 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar' |
|---|
| 1424 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1425 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1426 | CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:)) |
|---|
| 1427 | ENDIF |
|---|
| 1428 | |
|---|
| 1429 | |
|---|
| 1430 | |
|---|
| 1431 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar' |
|---|
| 1432 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1433 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1434 | CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:)) |
|---|
| 1435 | ENDIF |
|---|
| 1436 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar' |
|---|
| 1437 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1438 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1439 | CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:)) |
|---|
| 1440 | ENDIF |
|---|
| 1441 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar' |
|---|
| 1442 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1443 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1444 | CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:)) |
|---|
| 1445 | ENDIF |
|---|
| 1446 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar' |
|---|
| 1447 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1448 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1449 | CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:)) |
|---|
| 1450 | ENDIF |
|---|
| 1451 | |
|---|
| 1452 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar' |
|---|
| 1453 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1454 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1455 | CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:)) |
|---|
| 1456 | ENDIF |
|---|
| 1457 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar' |
|---|
| 1458 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1459 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1460 | CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:)) |
|---|
| 1461 | ENDIF |
|---|
| 1462 | |
|---|
| 1463 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar' |
|---|
| 1464 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1465 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1466 | CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:)) |
|---|
| 1467 | ENDIF |
|---|
| 1468 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar' |
|---|
| 1469 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1470 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1471 | CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:)) |
|---|
| 1472 | ENDIF |
|---|
| 1473 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar' |
|---|
| 1474 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1475 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1476 | CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:)) |
|---|
| 1477 | ENDIF |
|---|
| 1478 | |
|---|
| 1479 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar' |
|---|
| 1480 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1481 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1482 | CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:)) |
|---|
| 1483 | ENDIF |
|---|
| 1484 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar' |
|---|
| 1485 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1486 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1487 | CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:)) |
|---|
| 1488 | ENDIF |
|---|
| 1489 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar' |
|---|
| 1490 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1491 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1492 | CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:)) |
|---|
| 1493 | ENDIF |
|---|
| 1494 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar' |
|---|
| 1495 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1496 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1497 | CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:)) |
|---|
| 1498 | ENDIF |
|---|
| 1499 | |
|---|
| 1500 | |
|---|
| 1501 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar' |
|---|
| 1502 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1503 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1504 | CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:)) |
|---|
| 1505 | ENDIF |
|---|
| 1506 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar' |
|---|
| 1507 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1508 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1509 | CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:)) |
|---|
| 1510 | ENDIF |
|---|
| 1511 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar' |
|---|
| 1512 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1513 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1514 | CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:)) |
|---|
| 1515 | ENDIF |
|---|
| 1516 | |
|---|
| 1517 | tmp_u_amp_2d_mat(:,:) = 0. |
|---|
| 1518 | tmp_v_amp_2d_mat(:,:) = 0. |
|---|
| 1519 | tmp_u_phi_2d_mat(:,:) = 0. |
|---|
| 1520 | tmp_v_phi_2d_mat(:,:) = 0. |
|---|
| 1521 | |
|---|
| 1522 | a_u_2d_mat(:,:) = 0. |
|---|
| 1523 | b_u_2d_mat(:,:) = 0. |
|---|
| 1524 | a_v_2d_mat(:,:) = 0. |
|---|
| 1525 | b_v_2d_mat(:,:) = 0. |
|---|
| 1526 | |
|---|
| 1527 | qmax_2d_mat(:,:) = 0. |
|---|
| 1528 | qmin_2d_mat(:,:) = 0. |
|---|
| 1529 | |
|---|
| 1530 | ecc_2d_mat(:,:) = 0. |
|---|
| 1531 | thetamax_2d_mat(:,:) =0. |
|---|
| 1532 | thetamin_2d_mat(:,:) = 0. |
|---|
| 1533 | |
|---|
| 1534 | Qc_2d_mat(:,:) = 0. |
|---|
| 1535 | Qac_2d_mat(:,:) = 0. |
|---|
| 1536 | gc_2d_mat(:,:) = 0. |
|---|
| 1537 | gac_2d_mat(:,:) = 0. |
|---|
| 1538 | |
|---|
| 1539 | Phi_Ua_2d_mat(:,:) = 0. |
|---|
| 1540 | dir_Ua_2d_mat(:,:) = 0. |
|---|
| 1541 | polarity_2d_mat(:,:) = 0. |
|---|
| 1542 | |
|---|
| 1543 | |
|---|
| 1544 | END DO |
|---|
| 1545 | IF(lwp) WRITE(numout,*) "diaharm_fast: Finshed postprocessing 2d velocity tidal parameters" |
|---|
| 1546 | ENDIF |
|---|
| 1547 | |
|---|
| 1548 | CALL FLUSH(numout) |
|---|
| 1549 | |
|---|
| 1550 | IF (ln_ana_uv3d) THEN |
|---|
| 1551 | IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess baroclinic velocity tidal parameters" |
|---|
| 1552 | CALL FLUSH(numout) |
|---|
| 1553 | DO jh=1,nb_ana |
|---|
| 1554 | |
|---|
| 1555 | |
|---|
| 1556 | tmp_u_amp_3d_mat(:,:,:) = 0. |
|---|
| 1557 | tmp_v_amp_3d_mat(:,:,:) = 0. |
|---|
| 1558 | tmp_u_phi_3d_mat(:,:,:) = 0. |
|---|
| 1559 | tmp_v_phi_3d_mat(:,:,:) = 0. |
|---|
| 1560 | |
|---|
| 1561 | a_u_3d_mat(:,:,:) = 0. |
|---|
| 1562 | b_u_3d_mat(:,:,:) = 0. |
|---|
| 1563 | a_v_3d_mat(:,:,:) = 0. |
|---|
| 1564 | b_v_3d_mat(:,:,:) = 0. |
|---|
| 1565 | |
|---|
| 1566 | qmax_3d_mat(:,:,:) = 0. |
|---|
| 1567 | qmin_3d_mat(:,:,:) = 0. |
|---|
| 1568 | |
|---|
| 1569 | ecc_3d_mat(:,:,:) = 0. |
|---|
| 1570 | thetamax_3d_mat(:,:,:) =0. |
|---|
| 1571 | thetamin_3d_mat(:,:,:) = 0. |
|---|
| 1572 | |
|---|
| 1573 | Qc_3d_mat(:,:,:) = 0. |
|---|
| 1574 | Qac_3d_mat(:,:,:) = 0. |
|---|
| 1575 | gc_3d_mat(:,:,:) = 0. |
|---|
| 1576 | gac_3d_mat(:,:,:) = 0. |
|---|
| 1577 | |
|---|
| 1578 | Phi_Ua_3d_mat(:,:,:) = 0. |
|---|
| 1579 | dir_Ua_3d_mat(:,:,:) = 0. |
|---|
| 1580 | polarity_3d_mat(:,:,:) = 0. |
|---|
| 1581 | |
|---|
| 1582 | |
|---|
| 1583 | DO jk=1,jpkm1 |
|---|
| 1584 | DO jj = 1, nlcj !- 1 |
|---|
| 1585 | DO ji = 1, nlci ! - 1 |
|---|
| 1586 | |
|---|
| 1587 | IF ( ((umask(ji,jj,jk) + umask(ji-1,jj,jk)) > 0 ) .AND. ((vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) > 0 ) ) THEN |
|---|
| 1588 | |
|---|
| 1589 | if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then |
|---|
| 1590 | 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)) |
|---|
| 1591 | 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)))) |
|---|
| 1592 | else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then |
|---|
| 1593 | tmp_u_amp = (amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) |
|---|
| 1594 | tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) |
|---|
| 1595 | else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then |
|---|
| 1596 | tmp_u_amp = (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) |
|---|
| 1597 | tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)) |
|---|
| 1598 | else |
|---|
| 1599 | cycle |
|---|
| 1600 | end if |
|---|
| 1601 | |
|---|
| 1602 | |
|---|
| 1603 | if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then |
|---|
| 1604 | 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)) |
|---|
| 1605 | 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)))) |
|---|
| 1606 | else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then |
|---|
| 1607 | tmp_v_amp = (amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) |
|---|
| 1608 | tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) |
|---|
| 1609 | else if ( (vmask(ji,jj,jk) == 0) .AND. (vmask(ji,jj-1,jk) == 1)) then |
|---|
| 1610 | tmp_v_amp = (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)) |
|---|
| 1611 | tmp_v_phi = (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)) |
|---|
| 1612 | else |
|---|
| 1613 | cycle |
|---|
| 1614 | end if |
|---|
| 1615 | |
|---|
| 1616 | |
|---|
| 1617 | |
|---|
| 1618 | a_u = tmp_U_amp * cos(tmp_U_phi) |
|---|
| 1619 | b_u = tmp_U_amp * sin(tmp_U_phi) |
|---|
| 1620 | a_v = tmp_V_amp * cos(tmp_V_phi) |
|---|
| 1621 | b_v = tmp_V_amp * sin(tmp_V_phi) |
|---|
| 1622 | |
|---|
| 1623 | 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) ) ) ) |
|---|
| 1624 | delta = twodelta/2. |
|---|
| 1625 | |
|---|
| 1626 | !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)) ) |
|---|
| 1627 | 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)) |
|---|
| 1628 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1629 | |
|---|
| 1630 | !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)) ) |
|---|
| 1631 | alpha2 = sqrt( tmpreal ) |
|---|
| 1632 | if (alpha2 < 0) alpha2 = 0 !CYCLE |
|---|
| 1633 | alpha= sqrt( alpha2 ) |
|---|
| 1634 | |
|---|
| 1635 | |
|---|
| 1636 | !major and minor axis of the ellipse |
|---|
| 1637 | qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 ) |
|---|
| 1638 | |
|---|
| 1639 | tmpreal = (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 |
|---|
| 1640 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1641 | !qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 ) ! but always positive. |
|---|
| 1642 | qmin = sqrt( tmpreal ) ! but always positive. |
|---|
| 1643 | |
|---|
| 1644 | !eccentricity of ellipse |
|---|
| 1645 | tmpreal = (qmax + qmin) |
|---|
| 1646 | if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE |
|---|
| 1647 | !ecc = (qmax - qmin)/(qmax + qmin) |
|---|
| 1648 | ecc = (qmax - qmin)/(tmpreal) |
|---|
| 1649 | |
|---|
| 1650 | ! Angle of major and minor ellipse |
|---|
| 1651 | thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) ) |
|---|
| 1652 | thetamin = thetamax + rpi/2. |
|---|
| 1653 | |
|---|
| 1654 | |
|---|
| 1655 | |
|---|
| 1656 | ! Rotary current components: Pugh A3.10 |
|---|
| 1657 | ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors |
|---|
| 1658 | ! so Qc = clockwise = anticyclonic = negative |
|---|
| 1659 | ! and Qac = anticlockwise = cyclonic = negative |
|---|
| 1660 | |
|---|
| 1661 | tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) |
|---|
| 1662 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1663 | !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)) ) |
|---|
| 1664 | Qc = 0.5*sqrt( tmpreal ) |
|---|
| 1665 | |
|---|
| 1666 | tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) |
|---|
| 1667 | if (tmpreal < 0) tmpreal = 0 !CYCLE |
|---|
| 1668 | !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)) ) |
|---|
| 1669 | Qac = 0.5*sqrt( tmpreal ) |
|---|
| 1670 | |
|---|
| 1671 | |
|---|
| 1672 | gc = atan2( ( ( tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & |
|---|
| 1673 | & ( (tmp_U_amp*cos( tmp_U_phi )) - (tmp_V_amp*sin( tmp_V_phi )) ) ) |
|---|
| 1674 | gac = atan2( ( ( -tmp_U_amp*sin( tmp_U_phi ) ) + (tmp_V_amp*cos( tmp_V_phi) ) ) , & |
|---|
| 1675 | & ( (tmp_U_amp*cos( tmp_U_phi )) + (tmp_V_amp*sin( tmp_V_phi )) ) ) |
|---|
| 1676 | |
|---|
| 1677 | !Pugh A3.2 |
|---|
| 1678 | Phi_Ua = -0.5*(gac - gc) |
|---|
| 1679 | dir_Ua = 0.5*(gac + gc) ! positive from x axis |
|---|
| 1680 | |
|---|
| 1681 | tmpreal = qmax |
|---|
| 1682 | !if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE |
|---|
| 1683 | if (tmpreal == 0) then |
|---|
| 1684 | polarity = 0 |
|---|
| 1685 | else |
|---|
| 1686 | polarity = (Qac - Qc)/qmax |
|---|
| 1687 | endif |
|---|
| 1688 | |
|---|
| 1689 | |
|---|
| 1690 | tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp |
|---|
| 1691 | tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp |
|---|
| 1692 | tmp_u_phi_3d_mat(ji,jj,jk) = tmp_u_phi |
|---|
| 1693 | tmp_v_phi_3d_mat(ji,jj,jk) = tmp_v_phi |
|---|
| 1694 | |
|---|
| 1695 | |
|---|
| 1696 | a_u_3d_mat(ji,jj,jk) = a_u |
|---|
| 1697 | b_u_3d_mat(ji,jj,jk) = b_u |
|---|
| 1698 | a_v_3d_mat(ji,jj,jk) = a_v |
|---|
| 1699 | b_v_3d_mat(ji,jj,jk) = b_v |
|---|
| 1700 | |
|---|
| 1701 | qmax_3d_mat(ji,jj,jk) = qmax |
|---|
| 1702 | qmin_3d_mat(ji,jj,jk) = qmin |
|---|
| 1703 | |
|---|
| 1704 | ecc_3d_mat(ji,jj,jk) = ecc |
|---|
| 1705 | thetamax_3d_mat(ji,jj,jk) = thetamax |
|---|
| 1706 | thetamin_3d_mat(ji,jj,jk) = thetamin |
|---|
| 1707 | |
|---|
| 1708 | Qc_3d_mat(ji,jj,jk) = Qc |
|---|
| 1709 | Qac_3d_mat(ji,jj,jk) = Qac |
|---|
| 1710 | gc_3d_mat(ji,jj,jk) = gc |
|---|
| 1711 | gac_3d_mat(ji,jj,jk) = gac |
|---|
| 1712 | |
|---|
| 1713 | Phi_Ua_3d_mat(ji,jj,jk) = Phi_Ua |
|---|
| 1714 | dir_Ua_3d_mat(ji,jj,jk) = dir_Ua |
|---|
| 1715 | polarity_3d_mat(ji,jj,jk) = polarity |
|---|
| 1716 | |
|---|
| 1717 | ENDIF |
|---|
| 1718 | END DO !ji |
|---|
| 1719 | END DO !jj |
|---|
| 1720 | END DO !jk |
|---|
| 1721 | |
|---|
| 1722 | |
|---|
| 1723 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d' |
|---|
| 1724 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1725 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1726 | CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:)) |
|---|
| 1727 | ENDIF |
|---|
| 1728 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d' |
|---|
| 1729 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1730 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1731 | CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:)) |
|---|
| 1732 | ENDIF |
|---|
| 1733 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d' |
|---|
| 1734 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1735 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1736 | CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:)) |
|---|
| 1737 | ENDIF |
|---|
| 1738 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d' |
|---|
| 1739 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1740 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1741 | CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:)) |
|---|
| 1742 | ENDIF |
|---|
| 1743 | |
|---|
| 1744 | |
|---|
| 1745 | |
|---|
| 1746 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d' |
|---|
| 1747 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1748 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1749 | CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:)) |
|---|
| 1750 | ENDIF |
|---|
| 1751 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d' |
|---|
| 1752 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1753 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1754 | CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:)) |
|---|
| 1755 | ENDIF |
|---|
| 1756 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d' |
|---|
| 1757 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1758 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1759 | CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:)) |
|---|
| 1760 | ENDIF |
|---|
| 1761 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d' |
|---|
| 1762 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1763 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1764 | CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:)) |
|---|
| 1765 | ENDIF |
|---|
| 1766 | |
|---|
| 1767 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d' |
|---|
| 1768 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1769 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1770 | CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:)) |
|---|
| 1771 | ENDIF |
|---|
| 1772 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d' |
|---|
| 1773 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1774 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1775 | CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:)) |
|---|
| 1776 | ENDIF |
|---|
| 1777 | |
|---|
| 1778 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d' |
|---|
| 1779 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1780 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1781 | CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:)) |
|---|
| 1782 | ENDIF |
|---|
| 1783 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d' |
|---|
| 1784 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1785 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1786 | CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:)) |
|---|
| 1787 | ENDIF |
|---|
| 1788 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d' |
|---|
| 1789 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1790 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1791 | CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:)) |
|---|
| 1792 | ENDIF |
|---|
| 1793 | |
|---|
| 1794 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d' |
|---|
| 1795 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1796 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1797 | CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:)) |
|---|
| 1798 | ENDIF |
|---|
| 1799 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d' |
|---|
| 1800 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1801 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1802 | CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:)) |
|---|
| 1803 | ENDIF |
|---|
| 1804 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d' |
|---|
| 1805 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1806 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1807 | CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:)) |
|---|
| 1808 | ENDIF |
|---|
| 1809 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d' |
|---|
| 1810 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1811 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1812 | CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:)) |
|---|
| 1813 | ENDIF |
|---|
| 1814 | |
|---|
| 1815 | |
|---|
| 1816 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d' |
|---|
| 1817 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1818 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1819 | CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:)) |
|---|
| 1820 | ENDIF |
|---|
| 1821 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d' |
|---|
| 1822 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1823 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1824 | CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:)) |
|---|
| 1825 | ENDIF |
|---|
| 1826 | tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d' |
|---|
| 1827 | IF( iom_use(TRIM(tmp_name)) ) THEN |
|---|
| 1828 | IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name) |
|---|
| 1829 | CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:)) |
|---|
| 1830 | ENDIF |
|---|
| 1831 | |
|---|
| 1832 | tmp_u_amp_3d_mat(:,:,:) = 0. |
|---|
| 1833 | tmp_v_amp_3d_mat(:,:,:) = 0. |
|---|
| 1834 | tmp_u_phi_3d_mat(:,:,:) = 0. |
|---|
| 1835 | tmp_v_phi_3d_mat(:,:,:) = 0. |
|---|
| 1836 | |
|---|
| 1837 | a_u_3d_mat(:,:,:) = 0. |
|---|
| 1838 | b_u_3d_mat(:,:,:) = 0. |
|---|
| 1839 | a_v_3d_mat(:,:,:) = 0. |
|---|
| 1840 | b_v_3d_mat(:,:,:) = 0. |
|---|
| 1841 | |
|---|
| 1842 | qmax_3d_mat(:,:,:) = 0. |
|---|
| 1843 | qmin_3d_mat(:,:,:) = 0. |
|---|
| 1844 | |
|---|
| 1845 | ecc_3d_mat(:,:,:) = 0. |
|---|
| 1846 | thetamax_3d_mat(:,:,:) = 0. |
|---|
| 1847 | thetamin_3d_mat(:,:,:) = 0. |
|---|
| 1848 | |
|---|
| 1849 | Qc_3d_mat(:,:,:) = 0. |
|---|
| 1850 | Qac_3d_mat(:,:,:) = 0. |
|---|
| 1851 | gc_3d_mat(:,:,:) = 0. |
|---|
| 1852 | gac_3d_mat(:,:,:) = 0. |
|---|
| 1853 | |
|---|
| 1854 | Phi_Ua_3d_mat(:,:,:) = 0. |
|---|
| 1855 | dir_Ua_3d_mat(:,:,:) = 0. |
|---|
| 1856 | polarity_3d_mat(:,:,:) = 0. |
|---|
| 1857 | |
|---|
| 1858 | |
|---|
| 1859 | END DO !jh |
|---|
| 1860 | |
|---|
| 1861 | |
|---|
| 1862 | |
|---|
| 1863 | |
|---|
| 1864 | |
|---|
| 1865 | IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess 3d velocity tidal parameters" |
|---|
| 1866 | ENDIF |
|---|
| 1867 | |
|---|
| 1868 | CALL FLUSH(numout) |
|---|
| 1869 | ENDIF |
|---|
| 1870 | |
|---|
| 1871 | |
|---|
| 1872 | CALL FLUSH(numout) |
|---|
| 1873 | |
|---|
| 1874 | ! to output tidal parameters, u and v on t grid |
|---|
| 1875 | ! |
|---|
| 1876 | ! !== standard Cd ==! |
|---|
| 1877 | ! DO jj = 2, jpjm1 |
|---|
| 1878 | ! DO ji = 2, jpim1 |
|---|
| 1879 | ! imk = k_mk(ji,jj) ! ocean bottom level at t-points |
|---|
| 1880 | ! zut = un(ji,jj,imk) + un(ji-1,jj,imk) ! 2 x velocity at t-point |
|---|
| 1881 | ! zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk) |
|---|
| 1882 | ! ! ! here pCd0 = mask*boost * drag |
|---|
| 1883 | ! pCdU(ji,jj) = - pCd0(ji,jj) * SQRT( 0.25 * ( zut*zut + zvt*zvt ) + pke0 ) |
|---|
| 1884 | ! END DO |
|---|
| 1885 | ! END DO |
|---|
| 1886 | |
|---|
| 1887 | |
|---|
| 1888 | |
|---|
| 1889 | IF (ln_diaharm_postproc_vel) THEN |
|---|
| 1890 | IF (ln_ana_uvbar) THEN |
|---|
| 1891 | |
|---|
| 1892 | DEALLOCATE(amp_u2d, amp_v2d, phi_u2d, phi_v2d ) |
|---|
| 1893 | DEALLOCATE(TA_u_off_t_uvbar, TA_v_off_t_uvbar ) |
|---|
| 1894 | DEALLOCATE(TA_u_off_uvbar, TA_v_off_uvbar ) |
|---|
| 1895 | |
|---|
| 1896 | DEALLOCATE(tmp_u_amp_2d_mat, tmp_v_amp_2d_mat, tmp_u_phi_2d_mat, tmp_v_phi_2d_mat ) |
|---|
| 1897 | |
|---|
| 1898 | DEALLOCATE(a_u_2d_mat, b_u_2d_mat, a_v_2d_mat, b_v_2d_mat ) |
|---|
| 1899 | DEALLOCATE(qmax_2d_mat, qmin_2d_mat, ecc_2d_mat ) |
|---|
| 1900 | DEALLOCATE(thetamax_2d_mat, thetamin_2d_mat, Qc_2d_mat, Qac_2d_mat) |
|---|
| 1901 | DEALLOCATE(gc_2d_mat, gac_2d_mat, Phi_Ua_2d_mat, dir_Ua_2d_mat) |
|---|
| 1902 | DEALLOCATE(polarity_2d_mat ) |
|---|
| 1903 | |
|---|
| 1904 | ENDIF |
|---|
| 1905 | IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 2d velocity tidal parameters" |
|---|
| 1906 | |
|---|
| 1907 | IF (ln_ana_uv3d) THEN |
|---|
| 1908 | |
|---|
| 1909 | DEALLOCATE(amp_u3d, amp_v3d, phi_u3d, phi_v3d ) |
|---|
| 1910 | DEALLOCATE(TA_u_off_t_uv3d, TA_v_off_t_uv3d ) |
|---|
| 1911 | DEALLOCATE(TA_u_off_uv3d, TA_v_off_uv3d ) |
|---|
| 1912 | |
|---|
| 1913 | DEALLOCATE(tmp_u_amp_3d_mat, tmp_v_amp_3d_mat, tmp_u_phi_3d_mat, tmp_v_phi_3d_mat ) |
|---|
| 1914 | |
|---|
| 1915 | DEALLOCATE(a_u_3d_mat, b_u_3d_mat, a_v_3d_mat, b_v_3d_mat ) |
|---|
| 1916 | DEALLOCATE(qmax_3d_mat, qmin_3d_mat, ecc_3d_mat ) |
|---|
| 1917 | DEALLOCATE(thetamax_3d_mat, thetamin_3d_mat, Qc_3d_mat, Qac_3d_mat) |
|---|
| 1918 | DEALLOCATE(gc_3d_mat, gac_3d_mat, Phi_Ua_3d_mat, dir_Ua_3d_mat) |
|---|
| 1919 | DEALLOCATE(polarity_3d_mat ) |
|---|
| 1920 | |
|---|
| 1921 | ENDIF |
|---|
| 1922 | IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 3d velocity tidal parameters" |
|---|
| 1923 | |
|---|
| 1924 | ENDIF |
|---|
| 1925 | |
|---|
| 1926 | |
|---|
| 1927 | CALL FLUSH(numout) |
|---|
| 1928 | ! |
|---|
| 1929 | END SUBROUTINE harm_ana_out |
|---|
| 1930 | ! |
|---|
| 1931 | SUBROUTINE harm_rst_write(kt) |
|---|
| 1932 | !!---------------------------------------------------------------------- |
|---|
| 1933 | !! *** ROUTINE harm_ana_init *** |
|---|
| 1934 | !! |
|---|
| 1935 | !! ** Purpose : To write out cummulated Tidal Harmomnic data to file for |
|---|
| 1936 | !! restarting |
|---|
| 1937 | !! |
|---|
| 1938 | !! ** Method : restart files will be dated by default |
|---|
| 1939 | !! |
|---|
| 1940 | !! ** input : |
|---|
| 1941 | !! |
|---|
| 1942 | !! ** Action : ... |
|---|
| 1943 | !! |
|---|
| 1944 | !! history : |
|---|
| 1945 | !! 0.0 ! 01-16 (Enda O'Dea) Original code |
|---|
| 1946 | !! ASSUMES dated file for rose , can change later to be more generic |
|---|
| 1947 | !!---------------------------------------------------------------------- |
|---|
| 1948 | INTEGER, INTENT(in) :: kt ! ocean time-step |
|---|
| 1949 | !! |
|---|
| 1950 | INTEGER :: jh, j2d, j3d |
|---|
| 1951 | CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character |
|---|
| 1952 | CHARACTER(LEN=50) :: clname ! ocean output restart file name |
|---|
| 1953 | CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file |
|---|
| 1954 | CHARACTER(LEN=250) :: clfinal ! full name |
|---|
| 1955 | |
|---|
| 1956 | !restart file |
|---|
| 1957 | DO j2d=1,nvar_2d |
|---|
| 1958 | CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) ) |
|---|
| 1959 | DO jh=1,nb_ana |
|---|
| 1960 | CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) ) |
|---|
| 1961 | CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) ) |
|---|
| 1962 | ENDDO |
|---|
| 1963 | ENDDO |
|---|
| 1964 | |
|---|
| 1965 | DO j3d=1,nvar_3d |
|---|
| 1966 | !JT CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) ) |
|---|
| 1967 | CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) ) |
|---|
| 1968 | DO jh=1,nb_ana |
|---|
| 1969 | CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) ) |
|---|
| 1970 | CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) ) |
|---|
| 1971 | ENDDO |
|---|
| 1972 | ENDDO |
|---|
| 1973 | |
|---|
| 1974 | IF(lwp) THEN |
|---|
| 1975 | IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt |
|---|
| 1976 | ELSE ; WRITE(clkt, '(i8.8)') kt |
|---|
| 1977 | ENDIF |
|---|
| 1978 | clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin" |
|---|
| 1979 | clpath = TRIM(cn_ocerst_outdir) |
|---|
| 1980 | IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' |
|---|
| 1981 | IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname |
|---|
| 1982 | |
|---|
| 1983 | WRITE(clfinal,'(a)') trim(clpath)//trim(clname) |
|---|
| 1984 | OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" ) |
|---|
| 1985 | WRITE(66) cc |
|---|
| 1986 | WRITE(66) anau |
|---|
| 1987 | WRITE(66) anav |
|---|
| 1988 | WRITE(66) anaf |
|---|
| 1989 | WRITE(66) fjulday_startharm |
|---|
| 1990 | CLOSE(66) |
|---|
| 1991 | WRITE(numout,*) '----------------------------' |
|---|
| 1992 | WRITE(numout,*) ' harm_rst_write: DONE ' |
|---|
| 1993 | WRITE(numout,*) cc |
|---|
| 1994 | WRITE(numout,*) anaf |
|---|
| 1995 | WRITE(numout,*) anau |
|---|
| 1996 | WRITE(numout,*) anav |
|---|
| 1997 | WRITE(numout,*) fjulday_startharm |
|---|
| 1998 | WRITE(numout,*) '----------------------------' |
|---|
| 1999 | ENDIF |
|---|
| 2000 | |
|---|
| 2001 | END SUBROUTINE harm_rst_write |
|---|
| 2002 | |
|---|
| 2003 | SUBROUTINE harm_rst_read |
|---|
| 2004 | !!---------------------------------------------------------------------- |
|---|
| 2005 | !! *** ROUTINE harm_ana_init *** |
|---|
| 2006 | !! |
|---|
| 2007 | !! ** Purpose : To read in cummulated Tidal Harmomnic data to file for |
|---|
| 2008 | !! restarting |
|---|
| 2009 | !! |
|---|
| 2010 | !! ** Method : |
|---|
| 2011 | !! |
|---|
| 2012 | !! ** input : |
|---|
| 2013 | !! |
|---|
| 2014 | !! ** Action : ... |
|---|
| 2015 | !! |
|---|
| 2016 | !! history : |
|---|
| 2017 | !! 0.0 ! 01-16 (Enda O'Dea) Original code |
|---|
| 2018 | !! ASSUMES dated file for rose , can change later to be more generic |
|---|
| 2019 | !!---------------------------------------------------------------------- |
|---|
| 2020 | CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character |
|---|
| 2021 | CHARACTER(LEN=50) :: clname ! ocean output restart file name |
|---|
| 2022 | CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file |
|---|
| 2023 | CHARACTER(LEN=250) :: clfinal ! full name |
|---|
| 2024 | INTEGER :: jh, j2d, j3d |
|---|
| 2025 | |
|---|
| 2026 | IF( nit000 > 999999999 ) THEN ; WRITE(clkt, * ) nit000-1 |
|---|
| 2027 | ELSE ; WRITE(clkt, '(i8.8)') nit000-1 |
|---|
| 2028 | ENDIF |
|---|
| 2029 | clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin" |
|---|
| 2030 | clpath = TRIM(cn_ocerst_outdir) |
|---|
| 2031 | IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' |
|---|
| 2032 | |
|---|
| 2033 | IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname |
|---|
| 2034 | |
|---|
| 2035 | DO j2d=1,nvar_2d |
|---|
| 2036 | CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) ) |
|---|
| 2037 | IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) ) |
|---|
| 2038 | DO jh=1,nb_ana |
|---|
| 2039 | CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) ) |
|---|
| 2040 | CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) ) |
|---|
| 2041 | ENDDO |
|---|
| 2042 | ENDDO |
|---|
| 2043 | |
|---|
| 2044 | DO j3d=1,nvar_3d |
|---|
| 2045 | !JT CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) ) |
|---|
| 2046 | CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) ) |
|---|
| 2047 | IF(lwp) WRITE(numout,*) "3D", j3d, m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) ) |
|---|
| 2048 | |
|---|
| 2049 | DO jh=1,nb_ana |
|---|
| 2050 | CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) ) |
|---|
| 2051 | CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) ) |
|---|
| 2052 | ENDDO |
|---|
| 2053 | ENDDO |
|---|
| 2054 | |
|---|
| 2055 | WRITE(clfinal,'(a)') trim(clpath)//trim(clname) |
|---|
| 2056 | OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" ) |
|---|
| 2057 | READ(66) cc |
|---|
| 2058 | READ(66) anau |
|---|
| 2059 | READ(66) anav |
|---|
| 2060 | READ(66) anaf |
|---|
| 2061 | READ(66) fjulday_startharm |
|---|
| 2062 | CLOSE(66) |
|---|
| 2063 | |
|---|
| 2064 | IF(lwp) THEN |
|---|
| 2065 | WRITE(numout,*) '----------------------------' |
|---|
| 2066 | WRITE(numout,*) ' Checking anaf is correct' |
|---|
| 2067 | WRITE(numout,*) cc |
|---|
| 2068 | WRITE(numout,*) anaf |
|---|
| 2069 | WRITE(numout,*) fjulday_startharm |
|---|
| 2070 | WRITE(numout,*) '----------------------------' |
|---|
| 2071 | ENDIF |
|---|
| 2072 | |
|---|
| 2073 | END SUBROUTINE harm_rst_read |
|---|
| 2074 | |
|---|
| 2075 | !!====================================================================== |
|---|
| 2076 | |
|---|
| 2077 | END MODULE diaharm_fast |
|---|
