Changeset 12453


Ignore:
Timestamp:
2020-02-25T16:29:34+01:00 (8 months ago)
Author:
jcastill
Message:

First implementation of the branch - compiling after merge

Location:
NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE
Files:
2 added
14 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/BDY/bdyini.F90

    r11715 r12453  
    167167      ! Check and write out namelist parameters 
    168168      ! ----------------------------------------- 
    169       IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   & 
    170          &                               ' and general open boundary condition are not compatible' ) 
     169!      IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   & 
     170!         &                               ' and general open boundary condition are not compatible' ) 
    171171 
    172172      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DIA/diaharm.F90

    r11715 r12453  
    55   !!====================================================================== 
    66   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code 
     7   !!  
     8   !!   NB: 2017-12 : add 3D harmonic analysis of velocities  
     9   !!                 integration of Maria Luneva's development  
     10   !!   'key_3Ddiaharm 
    711   !!---------------------------------------------------------------------- 
    812   USE oce             ! ocean dynamics and tracers variables 
     
    1317   USE sbctide         ! Tidal forcing or not 
    1418   ! 
     19# if defined key_3Ddiaharm  
     20   USE zdf_oce  
     21#endif  
     22   ! 
    1523   USE in_out_manager  ! I/O units 
    1624   USE iom             ! I/0 library 
     
    3341   INTEGER         ::   nb_ana        ! Number of harmonics to analyse 
    3442 
    35    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
    36    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp 
    37    REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut   , vt   , ft 
    38    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta , out_u, out_v 
     43   INTEGER , ALLOCATABLE, DIMENSION(:)           ::   name 
     44   REAL(wp), ALLOCATABLE, DIMENSION(:)           ::   ana_freq, ut, vt, ft  
     45# if defined key_3Ddiaharm  
     46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)   ::   ana_temp  
     47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   out_eta, out_u, out_v, out_w, out_dzi  
     48# else  
     49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   ana_temp  
     50   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   out_eta, out_u, out_v  
     51# endif 
    3952 
    4053   INTEGER ::   ninco, nsparse 
     
    7689         WRITE(numout,*) 
    7790         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' 
     91# if defined key_3Ddiaharm  
     92         WRITE(numout,*) '  - 3D harmonic analysis of currents activated (key_3Ddiaharm)'  
     93#endif 
    7894         WRITE(numout,*) '~~~~~~~ ' 
    7995      ENDIF 
     
    155171         ! Initialize temporary arrays: 
    156172         ! ---------------------------- 
    157          ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    158          ana_temp(:,:,:,:) = 0._wp 
     173# if defined key_3Ddiaharm  
     174         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) )  
     175         ana_temp(:,:,:,:,:) = 0._wp  
     176# else  
     177         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3      ) )  
     178         ana_temp(:,:,:,:  ) = 0._wp  
     179#endif 
    159180 
    160181      ENDIF 
     
    175196      ! 
    176197      INTEGER  :: ji, jj, jh, jc, nhc 
     198# if defined key_3Ddiaharm  
     199      INTEGER  :: jk  
     200# endif 
    177201      REAL(wp) :: ztime, ztemp 
    178202      !!-------------------------------------------------------------------- 
     
    190214                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
    191215                  ! 
     216               ! ssh, ub, vb are stored at the last level of 5d array 
    192217               DO jj = 1,jpj 
    193218                  DO ji = 1,jpi 
    194                      ! Elevation 
    195                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)         
    196                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
    197                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
     219                     ! Elevation and currents  
     220# if defined key_3Ddiaharm  
     221                     ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj)          
     222                     ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj)  
     223                     ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)  
     224  
     225                     ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk)                               &  
     226                   &                              + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj)  
     227                     ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk)                               &   
     228                   &                              + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj)  
     229# else  
     230                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)          
     231                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)  
     232                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)  
     233# endif 
    198234                  END DO 
    199235               END DO 
    200236               ! 
     237# if defined key_3Ddiaharm  
     238! 3d velocity and density:  
     239             DO jk=1,jpk-1  
     240               DO jj = 1,jpj  
     241                  DO ji = 1,jpi  
     242                     ! density and velocity  
     243                     ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk)  
     244                     ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) &  
     245                &                                          *umask(ji,jj,jk)  
     246                     ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) &  
     247                &                                          *vmask(ji,jj,jk)   
     248                     ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk)  
     249   
     250                     ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1))/max(rn2(ji,jj,jk),1.e-8_wp)  
     251                  END DO  
     252               END DO  
     253             ENDDO  
     254# endif 
    201255            END DO 
    202256         END DO 
     
    221275      !!-------------------------------------------------------------------- 
    222276      INTEGER :: ji, jj, jh, jc, jn, nhan, jl 
     277# if defined key_3Ddiaharm  
     278      INTEGER :: jk  
     279# endif 
    223280      INTEGER :: ksp, kun, keq 
    224281      REAL(wp) :: ztime, ztime_ini, ztime_end 
     
    234291      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis 
    235292      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis 
     293 
     294# if defined key_3Ddiaharm  
     295      ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana),   &  
     296         &      out_u  (jpi,jpj,jpk,2*nb_ana),   &  
     297         &      out_v  (jpi,jpj,jpk,2*nb_ana),   &  
     298         &      out_w  (jpi,jpj,jpk,2*nb_ana),   &  
     299         &      out_dzi(jpi,jpj,jpk,2*nb_ana) )  
     300# else  
     301      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &  
     302         &      out_u  (jpi,jpj,2*nb_ana),   &  
     303         &      out_v  (jpi,jpj,2*nb_ana)  )  
     304# endif  
     305  
     306      IF(lwp) WRITE(numout,*) 'ANA F OLD', ft   
     307      IF(lwp) WRITE(numout,*) 'ANA U OLD', ut  
     308      IF(lwp) WRITE(numout,*) 'ANA V OLD', vt 
    236309 
    237310      ninco = 2*nb_ana 
     
    257330      nsparse = ksp 
    258331 
    259       ! Elevation: 
     332      ! Density and Elevation:  
     333# if defined key_3Ddiaharm  
     334    DO jk=1,jpk  
     335# endif 
    260336      DO jj = 1, jpj 
    261337         DO ji = 1, jpi 
     
    265341               DO jc = 1, 2 
    266342                  kun = kun + 1 
     343# if defined key_3Ddiaharm  
     344                  ztmp4(kun)=ana_temp(ji,jj,kun,1,jk)  
     345# else  
    267346                  ztmp4(kun)=ana_temp(ji,jj,kun,1) 
     347# endif 
    268348               END DO 
    269349            END DO 
     
    278358         END DO 
    279359      END DO 
    280  
    281       ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &  
    282          &      out_u  (jpi,jpj,2*nb_ana),   & 
    283          &      out_v  (jpi,jpj,2*nb_ana)  ) 
    284360 
    285361      DO jj = 1, jpj 
     
    288364               X1 = ana_amp(ji,jj,jh,1) 
    289365               X2 =-ana_amp(ji,jj,jh,2) 
    290                out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj) 
    291                out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj) 
    292             END DO 
    293          END DO 
    294       END DO 
    295  
    296       ! ubar: 
     366# if defined key_3Ddiaharm  
     367               out_eta(ji,jj,jk,jh       ) = X1 * tmask_i(ji,jj)  
     368               out_eta(ji,jj,jk,jh+nb_ana) = X2 * tmask_i(ji,jj)  
     369# else  
     370               out_eta(ji,jj   ,jh       ) = X1 * tmask_i(ji,jj)  
     371               out_eta(ji,jj   ,jh+nb_ana) = X2 * tmask_i(ji,jj)  
     372# endif  
     373            END DO  
     374         END DO  
     375      END DO  
     376  
     377      ! u-component of velocity 
    297378      DO jj = 1, jpj 
    298379         DO ji = 1, jpi 
     
    302383               DO jc = 1,2 
    303384                  kun = kun + 1 
     385# if defined key_3Ddiaharm  
     386                  ztmp4(kun)=ana_temp(ji,jj,kun,2,jk)  
     387# else 
    304388                  ztmp4(kun)=ana_temp(ji,jj,kun,2) 
     389# endif 
    305390               END DO 
    306391            END DO 
     
    322407               X1= ana_amp(ji,jj,jh,1) 
    323408               X2=-ana_amp(ji,jj,jh,2) 
    324                out_u(ji,jj,       jh) = X1 * ssumask(ji,jj) 
    325                out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 
     409# if defined key_3Ddiaharm  
     410               out_u(ji,jj,jk,       jh) = X1 * ssumask(ji,jj)  
     411               out_u(ji,jj,jk,nb_ana+jh) = X2 * ssumask(ji,jj)  
     412# else  
     413               out_u(ji,jj,          jh) = X1 * ssumask(ji,jj)  
     414               out_u(ji,jj,   nb_ana+jh) = X2 * ssumask(ji,jj)  
     415# endif 
    326416            ENDDO 
    327417         ENDDO 
    328418      ENDDO 
    329419 
    330       ! vbar: 
     420      ! v- velocity 
    331421      DO jj = 1, jpj 
    332422         DO ji = 1, jpi 
     
    336426               DO jc = 1,2 
    337427                  kun = kun + 1 
     428# if defined key_3Ddiaharm  
     429                  ztmp4(kun)=ana_temp(ji,jj,kun,3,jk)  
     430# else 
    338431                  ztmp4(kun)=ana_temp(ji,jj,kun,3) 
     432# endif 
    339433               END DO 
    340434            END DO 
     
    356450               X1=ana_amp(ji,jj,jh,1) 
    357451               X2=-ana_amp(ji,jj,jh,2) 
    358                out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj) 
    359                out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 
    360             END DO 
    361          END DO 
    362       END DO 
     452# if defined key_3Ddiaharm  
     453               out_v(ji,jj,jk,       jh)=X1 * ssvmask(ji,jj)  
     454               out_v(ji,jj,jk,nb_ana+jh)=X2 * ssvmask(ji,jj)  
     455# else  
     456               out_v(ji,jj,          jh)=X1 * ssvmask(ji,jj)  
     457               out_v(ji,jj,   nb_ana+jh)=X2 * ssvmask(ji,jj)  
     458# endif  
     459            END DO  
     460         END DO  
     461      END DO  
     462  
     463# if defined key_3Ddiaharm  
     464      ! w- velocity  
     465      DO jj = 1, jpj  
     466         DO ji = 1, jpi  
     467            ! Fill input array  
     468            kun=0  
     469            DO jh = 1,nb_ana  
     470               DO jc = 1,2  
     471                  kun = kun + 1  
     472                  ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)  
     473               END DO  
     474            END DO  
     475  
     476            CALL SUR_DETERMINE(jj+1)  
     477  
     478            ! Fill output array  
     479            DO jh = 1, nb_ana  
     480               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)  
     481               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)  
     482            END DO  
     483  
     484         END DO  
     485      END DO  
     486  
     487      DO jj = 1, jpj  
     488         DO ji = 1, jpi  
     489            DO jh = 1, nb_ana  
     490               X1=ana_amp(ji,jj,jh,1)  
     491               X2=-ana_amp(ji,jj,jh,2)  
     492               out_w(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)  
     493               out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)  
     494            END DO  
     495         END DO  
     496      END DO  
     497  
     498       ! dzi- isopycnal displacements  
     499      DO jj = 1, jpj  
     500         DO ji = 1, jpi  
     501            ! Fill input array  
     502            kun=0  
     503            DO jh = 1,nb_ana  
     504               DO jc = 1,2  
     505                  kun = kun + 1  
     506                  ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)  
     507               END DO  
     508            END DO  
     509  
     510            CALL SUR_DETERMINE(jj+1)  
     511  
     512            ! Fill output array  
     513            DO jh = 1, nb_ana  
     514               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)  
     515               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)  
     516            END DO  
     517  
     518         END DO  
     519      END DO  
     520  
     521      DO jj = 1, jpj  
     522         DO ji = 1, jpi  
     523            DO jh = 1, nb_ana  
     524               X1=ana_amp(ji,jj,jh,1)  
     525               X2=-ana_amp(ji,jj,jh,2)  
     526               out_dzi(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)  
     527               out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)  
     528            END DO  
     529         END DO  
     530      END DO  
     531  
     532   ENDDO ! jk  
     533# endif 
    363534      ! 
    364535      CALL dia_wri_harm ! Write results in files 
     
    379550         cdfile_name_V         ! name of the file created (V-points) 
    380551      INTEGER  ::   jh 
    381       !!---------------------------------------------------------------------- 
     552 
     553# if defined key_3Ddiaharm  
     554      CHARACTER(LEN=lc) :: cdfile_name_W         ! name of the file created (W-points)  
     555      INTEGER  :: jk  
     556      REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im   
     557      REAL(WP), ALLOCATABLE, DIMENSION (:,:)   :: z2real, z2im        
     558# endif  
     559!!----------------------------------------------------------------------  
     560  
     561#if defined key_dimgout  
     562      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'  
     563      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'  
     564      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'  
     565#   if defined key_3Ddiaharm  
     566      cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'  
     567#   endif  
     568#endif 
    382569 
    383570      IF(lwp) WRITE(numout,*) '  ' 
    384571      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
     572#if defined key_dimgout  
     573      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)  
     574      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)  
     575      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)  
     576#   if defined key_3Ddiaharm  
     577      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_W)  
     578#   endif  
     579#endif 
    385580      IF(lwp) WRITE(numout,*) '  ' 
    386581 
    387       ! A) Elevation 
     582# if defined key_3Ddiaharm  
     583      ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))  
     584# endif  
     585  
     586      ! A) density and elevation 
    388587      !///////////// 
    389588      ! 
     589#if defined key_dimgout  
     590      cltext='density amplitude and phase; elevation is level=jpk '  
     591      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')  
     592#else  
     593#   if defined key_3Ddiaharm  
     594      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
     595#   endif 
    390596      DO jh = 1, nb_ana 
     597#   if defined key_3Ddiaharm  
     598        DO jk=1,jpkm1  
     599          z3real(:,:,jk)=out_eta(:,:,jk,jh)  
     600          z3im  (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)  
     601        ENDDO  
     602      z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)  
     603      CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )  
     604      CALL iom_put( TRIM(tname(jh))//'y_ro', z3im  (:,:,:) )  
     605      CALL iom_put( TRIM(tname(jh))//'x'   , z2real(:,:  ) )  
     606      CALL iom_put( TRIM(tname(jh))//'y'   , z2im  (:,:  ) )  
     607#   else   
     608      WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh)) 
    391609      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 
    392610      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) ) 
    393       END DO 
    394  
    395       ! B) ubar 
     611#   endif  
     612      END DO  
     613#endif  
     614  
     615      ! B) u 
    396616      !///////// 
    397617      ! 
     618#if defined key_dimgout  
     619      cltext='3d u amplitude and phase; ubar is the last level'  
     620      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')  
     621#else  
     622#   if defined key_3Ddiaharm  
     623      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
     624#   endif 
    398625      DO jh = 1, nb_ana 
    399       CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 
    400       CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) ) 
    401       END DO 
    402  
    403       ! C) vbar 
     626#   if defined key_3Ddiaharm  
     627        DO jk=1,jpkm1  
     628          z3real(:,:,jk)=out_u(:,:,jk,jh)  
     629          z3im  (:,:,jk)=out_u(:,:,jk,jh+nb_ana)  
     630        ENDDO  
     631        z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)  
     632        CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )  
     633        CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:)  )  
     634        CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )  
     635        CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:)  )  
     636        z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)  
     637        CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )  
     638        CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:)  )  
     639#   else  
     640        CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )  
     641        CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )  
     642#   endif  
     643      END DO  
     644#endif  
     645  
     646      ! C) v 
    404647      !///////// 
    405648      ! 
     649#if defined key_dimgout  
     650      cltext='3d v amplitude and phase; vbar is the last level'  
     651      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')  
     652#else  
     653#   if defined key_3Ddiaharm  
     654      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp  
     655#   endif 
    406656      DO jh = 1, nb_ana 
    407          CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) ) 
    408          CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 
    409       END DO 
    410       ! 
     657#   if defined key_3Ddiaharm  
     658        DO jk=1,jpkm1  
     659          z3real(:,:,jk)=out_v(:,:,jk,jh)  
     660          z3im  (:,:,jk)=out_v(:,:,jk,jh+nb_ana)  
     661        ENDDO  
     662        z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)  
     663        CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )  
     664        CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:)  )  
     665        CALL iom_put( TRIM(tname(jh))//'x_v2d'  , z2real(:,:) )  
     666        CALL iom_put( TRIM(tname(jh))//'y_v2d'  , z2im (:,:)  )  
     667        z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)  
     668        CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )  
     669        CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:)  )  
     670#   else  
     671         CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )  
     672         CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )  
     673#   endif  
     674       END DO  
     675  
     676#endif  
     677      ! D) w  
     678# if defined key_3Ddiaharm  
     679#   if defined key_dimgout  
     680      cltext='3d w amplitude and phase; vort_baro is the last level'  
     681      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')  
     682#   else  
     683      DO jh = 1, nb_ana  
     684        DO jk=1,jpkm1  
     685         z3real(:,:,jk)=out_w(:,:,jk,jh)  
     686         z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)  
     687        ENDDO  
     688        CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )  
     689        CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )  
     690      END DO  
     691#   endif  
     692  
     693!       E) dzi + tau_bot  
     694#   if defined key_dimgout  
     695      cltext='dzi=g*ro/N2 amplitude and phase'  
     696      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')  
     697#   else  
     698      DO jh = 1, nb_ana  
     699        DO jk=1,jpkm1  
     700         z3real(:,:,jk)=out_dzi(:,:,jk,jh)  
     701         z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)  
     702        ENDDO  
     703        CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )  
     704        CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )  
     705      END DO  
     706#   endif  
     707# endif   
     708  
     709      !  
     710# if defined key_3Ddiaharm  
     711   DEALLOCATE(z3real, z3im, z2real,z2im)  
     712# endif 
     713 
    411714   END SUBROUTINE dia_wri_harm 
    412715 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DOM/dtatsd.F90

    r11715 r12453  
    2121   ! 
    2222   USE in_out_manager  ! I/O manager 
     23   USE iom             ! IOM library 
    2324   USE lib_mpp         ! MPP library 
    2425 
     
    3132   !                                  !!* namtsd  namelist : Temperature & Salinity Data * 
    3233   LOGICAL , PUBLIC ::   ln_tsd_init   !: T & S data flag 
     34   LOGICAL , PUBLIC ::   ln_tsd_interp !: vertical interpolation flag 
    3335   LOGICAL , PUBLIC ::   ln_tsd_dmp    !: internal damping toward input data flag 
    3436 
    3537   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     38   INTEGER                              ::   jpk_init, inum_dta  
     39   INTEGER                              ::   id, linum   ! local integers  
     40   INTEGER                              ::   zdim(4) 
    3641 
    3742   !!---------------------------------------------------------------------- 
     
    5358      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used 
    5459      ! 
    55       INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers 
     60      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers 
    5661      !! 
    5762      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files 
    58       TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read 
    59       TYPE(FLD_N)                   ::   sn_tem, sn_sal 
     63      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read  
     64      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk 
    6065      !! 
    61       NAMELIST/namtsd/   ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal 
     66      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_interp, ln_tsd_dmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk 
    6267      !!---------------------------------------------------------------------- 
    6368      ! 
    6469      !  Initialisation 
    65       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
     70      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
    6671      ! 
    6772      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
     
    8085         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    8186         WRITE(numout,*) '   Namelist namtsd' 
    82          WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init = ', ln_tsd_init 
    83          WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp  = ', ln_tsd_dmp 
     87         WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init   = ', ln_tsd_init 
     88         WRITE(numout,*) '      Interpolation of initial conditions in the vertical  ln_tsd_interp = ', ln_tsd_interp 
     89         WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp    = ', ln_tsd_dmp 
    8490         WRITE(numout,*) 
    8591         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN 
     
    94100         ln_tsd_init = .FALSE. 
    95101      ENDIF 
     102      IF( ln_tsd_interp .AND. ln_tsd_dmp ) THEN  
     103            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN  
     104      ENDIF  
     105      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN  
     106            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN  
     107      ENDIF 
    96108      ! 
    97109      !                             ! allocate the arrays (if necessary) 
    98110      IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 
    99111         ! 
    100          ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     112         IF( ln_tsd_interp ) THEN  
     113           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information  
     114         ELSE  
     115           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 )   
     116         ENDIF 
    101117         IF( ierr0 > 0 ) THEN 
    102118            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN 
    103119         ENDIF 
    104120         ! 
    105                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
    106          IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
    107                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
    108          IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    109          ! 
    110          IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 
     121         IF( ln_tsd_interp ) THEN  
     122            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )   
     123            id = iom_varid( inum_dta, sn_dep%clvar, zdim )  
     124            jpk_init = zdim(3)  
     125            IF(lwp) WRITE(numout,*) 'Dimension of vertical coordinate in ICs: ', jpk_init  
     126            CALL iom_close( inum_dta )   ! Close the input file  
     127            !  
     128                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )  
     129            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )  
     130                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )  
     131            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )    
     132                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )  
     133                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )  
     134         ELSE  
     135                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )  
     136            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )  
     137                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )  
     138            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )    
     139         ENDIF ! ln_tsd_interp  
     140  
     141         !  
     142         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN 
    111143            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN 
    112144         ENDIF 
    113145         !                         ! fill sf_tsd with sn_tem & sn_sal and control print 
    114146         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal 
     147         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk 
    115148         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print ) 
    116149         ! 
     
    138171      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data 
    139172      ! 
    140       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
    141       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers 
     173      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indices  
     174      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers 
    142175      REAL(wp)::   zl, zi                             ! local scalars 
    143       REAL(wp), DIMENSION(jpk) ::  ztp, zsp   ! 1D workspace 
    144176      !!---------------------------------------------------------------------- 
    145177      ! 
     
    176208!!gm end 
    177209      ! 
    178       ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask 
    179       ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  
    180       ! 
    181       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    182          ! 
    183          IF( kt == nit000 .AND. lwp )THEN 
    184             WRITE(numout,*) 
    185             WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh' 
    186          ENDIF 
    187          ! 
    188          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    189             DO ji = 1, jpi 
    190                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     210      IF( kt == nit000 .AND. lwp )THEN  
     211         WRITE(numout,*)  
     212         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'  
     213      ENDIF  
     214      !  
     215      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable  
     216      !  
     217         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points  
     218            DO jj= 1, jpj  
     219               DO ji= 1, jpi 
    191220                  zl = gdept_0(ji,jj,jk) 
    192                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    193                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    194                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    195                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    196                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    197                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    198                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    199                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    200                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    201                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    202                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    203                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     221                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data  
     222                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1)   
     223                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)  
     224                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data  
     225                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)  
     226                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)  
     227                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1  
     228                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)  
     229                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down  
     230                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)  
     231                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)  
     232                        ENDIF  
     233                        IF((zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN  
     234                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &  
     235                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))  
     236                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &  
     237                        &                          (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi  
     238                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &  
     239                        &                          (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi 
    204240                        ENDIF 
    205241                     END DO 
    206242                  ENDIF 
    207243               END DO 
    208                DO jk = 1, jpkm1 
    209                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    210                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    211                END DO 
    212                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    213                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
    214244            END DO 
    215245         END DO 
    216246         !  
     247         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)  
     248         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:) 
    217249      ELSE                                !==   z- or zps- coordinate   ==! 
    218250         !                              
    219          ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask 
    220          ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:) 
     251         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)  * tmask(:,:,:)  ! Mask  
     252         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:) 
    221253         ! 
    222254         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
     
    248280                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure 
    249281         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta ) 
     282         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure  
     283         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure  
    250284                                        DEALLOCATE( sf_tsd              )     ! the structure itself 
    251285      ENDIF 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90

    r11715 r12453  
    3333   USE dynadv         ! dynamics: vector invariant versus flux form 
    3434   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme 
     35   USE dynspg 
    3536   USE domvvl         ! variable volume 
    3637   USE bdy_oce   , ONLY: ln_bdy 
     
    178179            vn(:,:,jk) = va(:,:,jk) 
    179180         END DO 
     181         ! limit velocities  
     182         IF (ln_ulimit) THEN  
     183            call dyn_limit_velocity (kt)  
     184         ENDIF 
    180185         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    181186!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
     
    210215               END DO 
    211216            END DO 
     217            ! limit velocities  
     218            IF (ln_ulimit) THEN  
     219               call dyn_limit_velocity (kt)  
     220            ENDIF 
    212221            !                             ! ================! 
    213222         ELSE                             ! Variable volume ! 
     
    275284                  END DO 
    276285               END DO 
     286               ! limit velocities  
     287               IF (ln_ulimit) THEN  
     288                  call dyn_limit_velocity (kt)  
     289               ENDIF 
    277290               ! 
    278291            ELSE                          ! Asselin filter applied on thickness weighted velocity 
     
    302315                  END DO 
    303316               END DO 
     317               ! limit velocities  
     318               IF (ln_ulimit) THEN  
     319                  call dyn_limit_velocity (kt)  
     320               ENDIF 
    304321               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
    305322               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     
    377394   END SUBROUTINE dyn_nxt 
    378395 
     396   SUBROUTINE dyn_limit_velocity (kt)  
     397      ! limits maxming vlaues of un and vn by volume courant number  
     398      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index  
     399      !  
     400      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
     401      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm  
     402  
     403      ! limit fluxes  
     404      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number  
     405      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow  
     406  
     407      DO jk = 1, jpkm1  
     408         DO jj = 1, jpjm1  
     409            DO ji = 1, jpim1  
     410               ! U direction  
     411               zzu = un(ji,jj,jk)  
     412               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)  
     413               ! ips is 1 if flow is positive othersize zero  
     414               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
     415               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
     416               ! idev is 1 if divergent flow otherwise zero  
     417               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )  
     418               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )  
     419               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
     420               zzcn = zzcn * zcn  
     421               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &  
     422                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
     423               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &  
     424                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
     425               ! limit currents  
     426               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
     427               ! V  direction  
     428               zzu = vn(ji,jj,jk)  
     429               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)  
     430               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
     431               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
     432               ! idev is 1 if divergent flow otherwise zero  
     433               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )  
     434               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )  
     435               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
     436               zzcn = zzcn * zcn  
     437               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &  
     438                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
     439               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &  
     440                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
     441               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
     442            ENDDO  
     443         ENDDO  
     444      ENDDO  
     445  
     446    END SUBROUTINE dyn_limit_velocity 
     447 
    379448   !!========================================================================= 
    380449END MODULE dynnxt 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynspg.F90

    r11715 r12453  
    3939   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...  
    4040 
     41   LOGICAL, PUBLIC ::  ln_ulimit  
     42   REAL(wp), PUBLIC :: cn_ulimit,cnn_ulimit 
     43 
    4144   !                       ! Parameter to control the surface pressure gradient scheme 
    4245   INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting) 
     
    191194      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   & 
    192195      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   & 
    193       &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 
     196      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha, & 
     197      &                    ln_ulimit, cn_ulimit, cnn_ulimit 
    194198      !!---------------------------------------------------------------------- 
    195199      ! 
     
    213217         WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp 
    214218         WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts 
     219         WRITE(numout,*) '     Limit velocities                       ln_ulimit     = ', ln_ulimit  
     220         WRITE(numout,*) '     Limit velocities      max inverse Courant number     = ', cn_ulimit  
     221         WRITE(numout,*) '     Limit velocities   multiplier for divergant flow     = ', cnn_ulimit 
    215222      ENDIF 
    216223      !                          ! Control of surface pressure gradient scheme options 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/sbctide.F90

    r11715 r12453  
    1616   USE ioipsl         ! NetCDF IPSL library 
    1717   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     18   USE bdytides 
    1819 
    1920   IMPLICIT NONE 
     
    108109 
    109110      DO jk = 1, nb_harmo 
    110          zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk) 
     111         ! love number now provides in tide namelist  
     112         zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk) 
    111113         DO ji = 1, jpi 
    112114            DO jj = 1, jpj 
     
    119121               IF    ( Wave(ntide(jk))%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat ) 
    120122               ELSEIF( Wave(ntide(jk))%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2 
     123               ! Add tide potential for long period tides  
     124               ELSEIF( Wave(ntide(jk))%nutide == 0 )  THEN  ;  zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp) 
    121125               ELSE                                         ;  zcs = 0._wp 
    122126               ENDIF 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/tide_mod.F90

    r11715 r12453  
    1616   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules 
    1717 
     18# if defined key_FES14_tides  
     19   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 34   !: maximum number of harmonic  
     20# else 
    1821   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic 
     22# endif 
    1923 
    2024   TYPE, PUBLIC ::    tide 
     
    4145 
    4246   SUBROUTINE tide_init_Wave 
     47# if defined key_FES14_tides  
     48#     include "tide_FES14.h90"  
     49# else 
    4350#     include "tide.h90" 
     51# endif 
    4452   END SUBROUTINE tide_init_Wave 
    4553 
     
    331339         zf = zf * zf1 * zf1 
    332340         ! 
     341      CASE( 20 )                 !==  formule 20,  compound waves ( 78 x 78 x 78 x 78 )  
     342         zf1 = nodal_factort(78)  
     343         zf  = zf1 * zf1 * zf1 * zf1  
     344         ! 
    333345      CASE( 73 )                 !==  formule 73 
    334346         zs = sin(sh_I) 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/tideini.F90

    r11715 r12453  
    3333   INTEGER , PUBLIC ::   kt_tide         !: 
    3434   REAL(wp), PUBLIC ::   rdttideramp     !: 
     35   REAL(wp), PUBLIC ::   dn_love_number  !: 
    3536   REAL(wp), PUBLIC ::   rn_scal_load    !: 
    3637   CHARACTER(lc), PUBLIC ::   cn_tide_load   !:  
     
    5455      ! 
    5556      NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, & 
    56                   &     ln_tide_ramp, rn_scal_load, rdttideramp, clname 
     57                  &     ln_tide_ramp, rn_scal_load, rdttideramp, dn_love_number, clname 
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
     
    7879            WRITE(numout,*) '         Read load potential from file           ln_read_load = ', ln_read_load 
    7980            WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp = ', ln_tide_ramp 
     81            WRITE(numout,*) '                                              dn_love_number  = ', dn_love_number 
    8082            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load = ', rn_scal_load 
    8183            WRITE(numout,*) '         Duration (days) of ramp                 rdttideramp  = ', rdttideramp 
     
    99101      END DO 
    100102      !        
     103      IF (ln_tide .and.lwp) WRITE(numout,*) '                                     nb_harmo     = ', nb_harmo 
     104 
    101105      ! Ensure that tidal components have been set in namelist_cfg 
    102106      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/USR/usrdef_istate.F90

    r11715 r12453  
    6262         DO jj = 1, jpj 
    6363            DO ji = 1, jpi 
    64                pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   & 
    65                     &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             & 
    66                     &           + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) )            & 
    67                     &           - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.)                        & 
    68                     &           + 7.  * (1500. - pdept(ji,jj,jk) ) / 1500.)                        & 
    69                     &           * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * ptmask(ji,jj,jk) 
    70  
    71                pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) )  & 
    72                     &         * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2                  & 
    73                     &         + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000.                 & 
    74                     &         - 1.62 * TANH( (pdept(ji,jj,jk) - 60.  ) / 650. )                    & 
    75                     &         + 0.2  * TANH( (pdept(ji,jj,jk) - 35.  ) / 100. )                    & 
    76                     &         + 0.2  * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) )                  & 
    77                     &         * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * ptmask(ji,jj,jk) 
     64               pts(ji,jj,jk,jp_tem) = 20._wp * ptmask(ji,jj,jk)  
     65               pts(ji,jj,jk,jp_sal) = 36.25_wp * ptmask(ji,jj,jk) 
    7866            END DO 
    7967         END DO 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/USR/usrdef_sbc.F90

    r11715 r12453  
    33   !!                     ***  MODULE  usrdef_sbc  *** 
    44   !! 
    5    !!                     ===  GYRE configuration  === 
     5   !!                  ===  WAD_TEST_CASES configuration  === 
    66   !! 
    77   !! User defined :   surface forcing of a user configuration 
     
    1111 
    1212   !!---------------------------------------------------------------------- 
    13    !!   usrdef_sbc    : user defined surface bounday conditions in GYRE case 
     13   !!   usrdef_sbc    : user defined surface bounday conditions in WAD_TEST_CASES case 
    1414   !!---------------------------------------------------------------------- 
    1515   USE oce            ! ocean dynamics and tracers 
     
    2121   USE lib_mpp        ! distribued memory computing library 
    2222   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    23    USE lib_fortran    ! 
     23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2424 
    2525   IMPLICIT NONE 
     
    4343      !!                    ***  ROUTINE usrdef_sbc  *** 
    4444      !!               
    45       !! ** Purpose :   provide at each time-step the GYRE surface boundary 
     45      !! ** Purpose :   provide at each time-step the surface boundary 
    4646      !!              condition, i.e. the momentum, heat and freshwater fluxes. 
    4747      !! 
    48       !! ** Method  :   analytical seasonal cycle for GYRE configuration. 
     48      !! ** Method  :   all 0 fields, for WAD_TEST_CASES case 
    4949      !!                CAUTION : never mask the surface stress field ! 
    5050      !! 
    51       !! ** Action  : - set the ocean surface boundary condition, i.e.    
     51      !! ** Action  : - set to ZERO all the ocean surface boundary condition, i.e. 
    5252      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx 
    5353      !! 
    54       !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000. 
    5554      !!---------------------------------------------------------------------- 
    5655      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    57       !! 
    58       INTEGER  ::   ji, jj                 ! dummy loop indices 
    59       INTEGER  ::   zyear0                 ! initial year  
    60       INTEGER  ::   zmonth0                ! initial month 
    61       INTEGER  ::   zday0                  ! initial day 
    62       INTEGER  ::   zday_year0             ! initial day since january 1st 
    63       REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle 
    64       REAL(wp) ::   ztime                  ! time in hour 
    65       REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january 
    66       REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january 
    67       REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january 
    68       REAL(wp) ::   ztaun                  ! intensity 
    69       REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar 
    70       REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star 
    71       REAL(wp) ::   zsumemp, zsurf 
    72       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3 
    73       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
    74       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
    75       REAL(wp) ::   zyydd                 ! number of days in one year 
    7656      !!--------------------------------------------------------------------- 
    77       zyydd = REAL(nyear_len(1),wp) 
    78  
    79       ! ---------------------------- ! 
    80       !  heat and freshwater fluxes  ! 
    81       ! ---------------------------- ! 
    82       !same temperature, E-P as in HAZELEGER 2000 
    83  
    84       zyear0     =   ndate0 / 10000                             ! initial year 
    85       zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month 
    86       zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30 
    87       zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360 
    88  
    89       ! current day (in hours) since january the 1st of the current year 
    90       ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours) 
    91          &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours) 
    92  
    93       ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours 
    94       ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours 
    95       ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours 
    96       ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours 
    97       !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours 
    98  
    99       ! amplitudes 
    100       zemp_S    = 0.7       ! intensity of COS in the South 
    101       zemp_N    = 0.8       ! intensity of COS in the North 
    102       zemp_sais = 0.1 
    103       zTstar    = 28.3      ! intemsity from 28.3 a -5 deg 
    104  
    105       ! 1/2 period between 21th June and 21th December and between 21th July and 21th January 
    106       zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi )  
    107       zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi ) 
    108  
    109       ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K) 
    110       zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s 
    111       DO jj = 1, jpj 
    112          DO ji = 1, jpi 
    113             ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3 
    114             ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period : 
    115             ! 64.5 in summer, 42.5 in winter 
    116             t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                & 
    117                &                    * COS( rpi * (gphit(ji,jj) - 5.)               & 
    118                &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) ) 
    119             ! 23.5 deg : tropics 
    120             qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) ) 
    121             qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj) 
    122             IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg 
    123                emp  (ji,jj) =   zemp_S * zconv   & 
    124                   &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & 
    125                   &         * ( 1 - zemp_sais / zemp_S * zcos_sais1) 
    126             ELSE 
    127                emp (ji,jj) =  - zemp_N * zconv   & 
    128                   &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & 
    129                   &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 ) 
    130             ENDIF 
    131          END DO 
    132       END DO 
    133  
    134       zsumemp = GLOB_SUM( 'usrdef_sbc', emp  (:,:)   )  
    135       zsurf   = GLOB_SUM( 'usrdef_sbc', tmask(:,:,1) )  
    136       zsumemp = zsumemp / zsurf         ! Default GYRE configuration 
    137  
    138       ! freshwater (mass flux) and update of qns with heat content of emp 
    139       emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average) 
    140       sfx (:,:) = 0.0_wp                                   ! no salt flux 
    141       qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST 
    142  
    143  
    144       ! ---------------------------- ! 
    145       !       momentum fluxes        ! 
    146       ! ---------------------------- ! 
    147       ! same wind as in Wico 
    148       !test date0 : ndate0 = 010203 
    149       zyear0  =   ndate0 / 10000 
    150       zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100 
    151       zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100 
    152       !Calculates nday_year, day since january 1st 
    153       zday_year0 = (zmonth0-1)*30.+zday0 
    154  
    155       !accumulates days of previous months of this year 
    156       ! day (in hours) since january the 1st 
    157       ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour 
    158          &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years 
    159       ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours 
    160       ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours 
    161       !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours 
    162  
    163       ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle 
    164       ztau = 0.105 / SQRT( 2. ) 
    165       ! seasonal oscillation intensity 
    166       ztau_sais = 0.015 
    167       ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi ) 
    168       DO jj = 1, jpj 
    169          DO ji = 1, jpi 
    170            ! domain from 15deg to 50deg and 1/2 period along 14deg 
    171            ! so 5/4 of half period with seasonal cycle 
    172            utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) ) 
    173            vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) ) 
    174          END DO 
    175       END DO 
    176  
    177       ! module of wind stress and wind speed at T-point 
    178       zcoef = 1. / ( zrhoa * zcdrag )  
    179       DO jj = 2, jpjm1 
    180          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    181             ztx = utau(ji-1,jj  ) + utau(ji,jj)  
    182             zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
    183             zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    184             taum(ji,jj) = zmod 
    185             wndm(ji,jj) = SQRT( zmod * zcoef ) 
    186          END DO 
    187       END DO 
    188       CALL lbc_lnk_multi( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. ) 
    189  
    190       ! ---------------------------------- ! 
    191       !  control print at first time-step  ! 
    192       ! ---------------------------------- ! 
    193       IF( kt == nit000 .AND. lwp ) THEN  
    194          WRITE(numout,*) 
    195          WRITE(numout,*)'usrdef_sbc_oce : analytical surface fluxes for GYRE configuration'                
    196          WRITE(numout,*)'~~~~~~~~~~~ '  
    197          WRITE(numout,*)'           nyear      = ', nyear 
    198          WRITE(numout,*)'           nmonth     = ', nmonth 
    199          WRITE(numout,*)'           nday       = ', nday 
    200          WRITE(numout,*)'           nday_year  = ', nday_year 
    201          WRITE(numout,*)'           ztime      = ', ztime 
    202          WRITE(numout,*)'           ztimemax   = ', ztimemax 
    203          WRITE(numout,*)'           ztimemin   = ', ztimemin 
    204          WRITE(numout,*)'           ztimemax1  = ', ztimemax1 
    205          WRITE(numout,*)'           ztimemin1  = ', ztimemin1 
    206          WRITE(numout,*)'           ztimemax2  = ', ztimemax2 
    207          WRITE(numout,*)'           ztimemin2  = ', ztimemin2 
    208          WRITE(numout,*)'           zyear0     = ', zyear0 
    209          WRITE(numout,*)'           zmonth0    = ', zmonth0 
    210          WRITE(numout,*)'           zday0      = ', zday0 
    211          WRITE(numout,*)'           zday_year0 = ', zday_year0 
    212          WRITE(numout,*)'           zyydd      = ', zyydd 
    213          WRITE(numout,*)'           zemp_S     = ', zemp_S 
    214          WRITE(numout,*)'           zemp_N     = ', zemp_N 
    215          WRITE(numout,*)'           zemp_sais  = ', zemp_sais 
    216          WRITE(numout,*)'           zTstar     = ', zTstar 
    217          WRITE(numout,*)'           zsumemp    = ', zsumemp 
    218          WRITE(numout,*)'           zsurf      = ', zsurf 
    219          WRITE(numout,*)'           ztrp       = ', ztrp 
    220          WRITE(numout,*)'           zconv      = ', zconv 
    221          WRITE(numout,*)'           ndastp     = ', ndastp 
    222          WRITE(numout,*)'           adatrj     = ', adatrj 
     57      !  
     58      IF( kt == nit000 ) THEN  
     59         !  
     60         IF(lwp) WRITE(numout,*)' usr_sbc : WAD_TEST_CASES case: NO surface forcing'  
     61         IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0'  
     62         !  
     63         utau(:,:) = 0._wp  
     64         vtau(:,:) = 0._wp  
     65         taum(:,:) = 0._wp  
     66         wndm(:,:) = 0._wp  
     67         !  
     68         emp (:,:) = 0._wp  
     69         sfx (:,:) = 0._wp  
     70         qns (:,:) = 0._wp  
     71         qsr (:,:) = 0._wp  
     72         ! 
    22373      ENDIF 
    22474      ! 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/par_oce.F90

    r11715 r12453  
    6464   INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature 
    6565   INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity 
     66   INTEGER, PUBLIC, PARAMETER ::   jp_dep = 3    !: indice for depth  
     67   INTEGER, PUBLIC, PARAMETER ::   jp_msk = 4    !: indice for mask 
    6668 
    6769   !!---------------------------------------------------------------------- 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/step.F90

    r11715 r12453  
    209209                         CALL dia_ar5 ( kstp )        ! ar5 diag 
    210210      IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     211      IF( lk_diaharm_fast )                           &  
     212            &            CALL dia_harm_fast( kstp )   ! Tidal harmonic analysis - restart and faster version 
    211213                         CALL dia_wri ( kstp )        ! ocean model: outputs 
    212214      ! 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/step_oce.F90

    r11715 r12453  
    8080   USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine) 
    8181   USE diaharm 
     82   USE diaharm_fast     ! harmonic analysis of tides (harm_ana routine) 
    8283   USE diacfl 
    8384   USE diaobs          ! Observation operator 
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/stpctl.F90

    r11715 r12453  
    115115         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max 
    116116      ENDIF 
    117       zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only) 
     117      zmax(2) = MAXVAL(  un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)  )         ! velocity max 
    118118      zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
    119119      zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
     
    149149      IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 
    150150         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
    151          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     151         &  zmax(2) >  100._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    152152         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
    153153         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
     
    156156         IF( lk_mpp .AND. ln_ctl ) THEN 
    157157            CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih  ) 
    158             CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  ) 
     158            CALL mpp_maxloc( 'stpctl', un*un + vn*vn    , umask (:,:,:), zzz, iu  ) 
    159159            CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 ) 
    160160            CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 ) 
    161161         ELSE 
    162162            ih(:)  = MAXLOC( ABS( sshn(:,:)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
    163             iu(:)  = MAXLOC( ABS( un  (:,:,:) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     163            iu(:)  = MAXLOC( un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)       ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    164164            is1(:) = MINLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    165165            is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     
    187187      ! 
    1881889100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5) 
    189 9200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
     1899200  FORMAT (' kt=',i8,'   Vel2  max: ',1pg11.4,', at  i j k: ',3i5) 
    1901909300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5) 
    1911919400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5) 
    192 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
     1929500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' Vel2_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
    193193      ! 
    194194   END SUBROUTINE stp_ctl 
Note: See TracChangeset for help on using the changeset viewer.