Changeset 5777


Ignore:
Timestamp:
2015-10-06T13:40:42+02:00 (5 years ago)
Author:
gm
Message:

#1593: LDF-ADV, III. Phasing of the improvements/simplifications of ADV & LDF momentum trends (see wiki page)

Location:
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO
Files:
8 deleted
93 edited
2 moved

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90

    r4624 r5777  
    7171         CALL fld_read( kt, nn_fsbc, sf_icedmp ) 
    7272         ! 
    73 !CDIR COLLAPSE 
    7473         hicif(:,:) = MAX( 0._wp,                     &        ! h >= 0         avoid spurious out of physical range 
    7574            &         hicif(:,:) - rdt_ice * resto_ice(:,:,1) * ( hicif(:,:) - sf_icedmp(jp_hicif)%fnow(:,:,1) )  )  
    76 !CDIR COLLAPSE 
    7775         frld (:,:) = MAX( 0._wp, MIN( 1._wp,         &        ! 0<= frld<=1    values which blow the run up 
    7876            &         frld (:,:) - rdt_ice * resto_ice(:,:,1) * ( frld (:,:) - sf_icedmp(jp_frld )%fnow(:,:,1) )  )  ) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90

    r5123 r5777  
    160160      !------------------------------------------------------------------- 
    161161 
    162 !CDIR NOVERRCHK 
    163162      DO jj = k_j1 , k_jpj-1 
    164 !CDIR NOVERRCHK 
    165163         DO ji = 1 , jpi 
    166164            ! only the sinus changes its sign with the hemisphere 
     
    245243         ! Computation of free drift field for free slip boundary conditions. 
    246244 
    247 !CDIR NOVERRCHK 
    248245         DO jj = k_j1, k_jpj-1 
    249 !CDIR NOVERRCHK 
    250246            DO ji = 1, fs_jpim1 
    251247               !- Rate of strain tensor. 
     
    401397iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
    402398            !                                                   ! ================ ! 
    403 !CDIR NOVERRCHK 
    404399            DO jj = k_j1+1, k_jpj-1 
    405 !CDIR NOVERRCHK 
    406400               DO ji = 2, fs_jpim1   ! NO vector opt. 
    407401                  ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90

    r5407 r5777  
    319319         ! 
    320320         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step) 
    321 !CDIR NOVERRCHK 
     321            ! 
    322322            DO jj = 1, jpj                               !* modulus of ice-ocean relative velocity at I-point 
    323 !CDIR NOVERRCHK 
    324323               DO ji = 1, jpi 
    325324                  zu_i  = u_ice(ji,jj) - u_oce(ji,jj)                   ! ice-ocean relative velocity at I-point 
     
    328327               END DO 
    329328            END DO 
    330 !CDIR NOVERRCHK 
    331329            DO jj = 1, jpjm1                             !* update the modulus of stress at ocean surface (T-point) 
    332 !CDIR NOVERRCHK 
    333330               DO ji = 1, jpim1   ! NO vector opt. 
    334331                  !                                               ! modulus of U_ice-U_oce at T-point 
     
    383380         ! 
    384381         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step) 
    385 !CDIR NOVERRCHK 
     382            ! 
    386383            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity at T-point 
    387 !CDIR NOVERRCHK 
    388384               DO ji = fs_2, fs_jpim1 
    389385                  zu_t  = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   ! 2*(U_ice-U_oce) at T-point 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90

    r5407 r5777  
    196196      !-------------------------------------------------------------------------- 
    197197 
    198       !CDIR NOVERRCHK 
    199198      DO jj = 1, jpj 
    200          !CDIR NOVERRCHK 
    201199         DO ji = 1, jpi 
    202200            zthsnice       = hsnif(ji,jj) + hicif(ji,jj) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r3625 r5777  
    134134      !--------------------------------------------------------------------- 
    135135       
    136 !CDIR NOVERRCHK 
    137136      DO ji = kideb , kiut 
    138137         iicefr       = 1 - MAX( 0, INT( SIGN( 1.5 * zone , zfrl_old(ji) - 1.0 + epsi13 ) ) ) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5737 r5777  
    1414 
    1515   !!---------------------------------------------------------------------- 
    16    !!   'key_asminc'   : Switch on the assimilation increment interface 
    17    !!---------------------------------------------------------------------- 
    1816   !!   asm_inc_init   : Initialize the increment arrays and IAU weights 
    1917   !!   calc_date      : Compute the calendar date YYYYMMDD on a given step 
     
    2826   USE domvvl           ! domain: variable volume level 
    2927   USE oce              ! Dynamics and active tracers defined in memory 
    30    USE ldfdyn_oce       ! ocean dynamics: lateral physics 
     28   USE ldfdyn           ! lateral diffusion: eddy viscosity coefficients 
    3129   USE eosbn2           ! Equation of state - in situ and potential density 
    3230   USE zpshde           ! Partial step : Horizontal Derivative 
     
    5654    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments 
    5755#endif 
    58    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
    59    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    60    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
    61    LOGICAL, PUBLIC :: ln_trainc = .FALSE.      !: No tracer (T and S) assimilation increments 
    62    LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments 
    63    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    64    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment 
    65    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
     56   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields 
     57   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment 
     58   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization 
     59   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments 
     60   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments 
     61   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment 
     62   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment 
     63   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check 
    6664   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
    67    INTEGER, PUBLIC :: nn_divdmp                !: Apply divergence damping filter nn_divdmp times 
     65   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times 
    6866 
    6967   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity 
     
    9088   !! * Substitutions 
    9189#  include "domzgr_substitute.h90" 
    92 #  include "ldfdyn_substitute.h90" 
    9390#  include "vectopt_loop_substitute.h90" 
    9491   !!---------------------------------------------------------------------- 
     
    139136      ! Read Namelist nam_asminc : assimilation increment interface 
    140137      !----------------------------------------------------------------------- 
    141       ln_seaiceinc = .FALSE. 
     138      ln_seaiceinc   = .FALSE. 
    142139      ln_temnofreeze = .FALSE. 
    143140 
     
    428425 
    429426      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 
    430  
    431          CALL wrk_alloc(jpi,jpj,hdiv)  
    432  
    433          DO  jt = 1, nn_divdmp 
    434  
     427         ! 
     428         CALL wrk_alloc( jpi,jpj,   hdiv )  
     429         ! 
     430         DO jt = 1, nn_divdmp 
     431            ! 
    435432            DO jk = 1, jpkm1 
    436  
    437433               hdiv(:,:) = 0._wp 
    438  
    439434               DO jj = 2, jpjm1 
    440435                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    447442                  END DO 
    448443               END DO 
    449  
    450444               CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change) 
    451  
     445               ! 
    452446               DO jj = 2, jpjm1 
    453447                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    454448                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji+1,jj) * hdiv(ji+1,jj)   & 
    455                                                                         - e1e2t(ji  ,jj) * hdiv(ji  ,jj) ) & 
    456                                                                       * r1_e1u(ji,jj) * umask(ji,jj,jk)  
     449                        &                                               - e1e2t(ji  ,jj) * hdiv(ji  ,jj) ) & 
     450                        &                                             * r1_e1u(ji,jj) * umask(ji,jj,jk)  
    457451                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1e2t(ji,jj+1) * hdiv(ji,jj+1)   & 
    458                                                                         - e1e2t(ji,jj  ) * hdiv(ji,jj  ) ) & 
    459                                                                       * r1_e2v(ji,jj) * vmask(ji,jj,jk)  
     452                        &                                               - e1e2t(ji,jj  ) * hdiv(ji,jj  ) ) & 
     453                        &                                             * r1_e2v(ji,jj) * vmask(ji,jj,jk)  
    460454                  END DO 
    461455               END DO 
    462  
    463456            END DO 
    464  
     457            ! 
    465458         END DO 
    466  
    467          CALL wrk_dealloc(jpi,jpj,hdiv)  
    468  
     459         ! 
     460         CALL wrk_dealloc( jpi,jpj,   hdiv )  
     461         ! 
    469462      ENDIF 
    470  
    471  
    472463 
    473464      !----------------------------------------------------------------------- 
     
    476467 
    477468      IF ( ln_asmdin ) THEN 
    478  
     469         ! 
    479470         ALLOCATE( t_bkg(jpi,jpj,jpk) ) 
    480471         ALLOCATE( s_bkg(jpi,jpj,jpk) ) 
     
    482473         ALLOCATE( v_bkg(jpi,jpj,jpk) ) 
    483474         ALLOCATE( ssh_bkg(jpi,jpj)   ) 
    484  
    485          t_bkg(:,:,:) = 0.0 
    486          s_bkg(:,:,:) = 0.0 
    487          u_bkg(:,:,:) = 0.0 
    488          v_bkg(:,:,:) = 0.0 
    489          ssh_bkg(:,:) = 0.0 
    490  
     475         ! 
     476         t_bkg(:,:,:) = 0._wp 
     477         s_bkg(:,:,:) = 0._wp 
     478         u_bkg(:,:,:) = 0._wp 
     479         v_bkg(:,:,:) = 0._wp 
     480         ssh_bkg(:,:) = 0._wp 
     481         ! 
    491482         !-------------------------------------------------------------------- 
    492483         ! Read from file the background state at analysis time 
    493484         !-------------------------------------------------------------------- 
    494  
     485         ! 
    495486         CALL iom_open( c_asmdin, inum ) 
    496  
     487         ! 
    497488         CALL iom_get( inum, 'rdastp', zdate_bkg )  
    498          
     489         ! 
    499490         IF(lwp) THEN 
    500491            WRITE(numout,*)  
    501             WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', & 
    502                &  NINT( zdate_bkg ) 
     492            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', NINT( zdate_bkg ) 
    503493            WRITE(numout,*) '~~~~~~~~~~~~' 
    504494         ENDIF 
    505  
     495         ! 
    506496         IF ( NINT( zdate_bkg ) /= iitdin_date ) & 
    507497            & CALL ctl_warn( ' Validity time of assimilation background state does', & 
    508498            &                ' not agree with Direct Initialization time' ) 
    509  
     499         ! 
    510500         IF ( ln_trainc ) THEN    
    511501            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg ) 
     
    514504            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:) 
    515505         ENDIF 
    516  
     506         ! 
    517507         IF ( ln_dyninc ) THEN    
    518508            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg ) 
     
    521511            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:) 
    522512         ENDIF 
    523          
     513         ! 
    524514         IF ( ln_sshinc ) THEN 
    525515            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg ) 
    526516            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1) 
    527517         ENDIF 
    528  
     518         ! 
    529519         CALL iom_close( inum ) 
    530  
     520         ! 
    531521      ENDIF 
    532522      ! 
     
    574564      ! If kt = kit000 - 1 then set the date to the restart date 
    575565      IF ( kt == kit000 - 1 ) THEN 
    576  
    577566         kdate = ndastp 
    578567         RETURN 
    579  
    580568      ENDIF 
    581569 
     
    646634      !! ** Action  :  
    647635      !!---------------------------------------------------------------------- 
    648       INTEGER, INTENT(IN) :: kt               ! Current time step 
    649       ! 
    650       INTEGER :: ji,jj,jk 
    651       INTEGER :: it 
     636      INTEGER, INTENT(IN) ::   kt   ! Current time step 
     637      ! 
     638      INTEGER  :: ji, jj, jk 
     639      INTEGER  :: it 
    652640      REAL(wp) :: zincwgt  ! IAU weight for current time step 
    653641      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 
    654642      !!---------------------------------------------------------------------- 
    655  
     643      ! 
    656644      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)  
    657645      ! used to prevent the applied increments taking the temperature below the local freezing point  
    658  
    659646      DO jk = 1, jpkm1 
    660647        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) ) 
    661648      END DO 
    662  
    663       IF ( ln_asmiau ) THEN 
    664  
    665          !-------------------------------------------------------------------- 
    666          ! Incremental Analysis Updating 
    667          !-------------------------------------------------------------------- 
    668  
     649         ! 
     650         !                             !-------------------------------------- 
     651      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating 
     652         !                             !-------------------------------------- 
     653         ! 
    669654         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    670  
     655            ! 
    671656            it = kt - nit000 + 1 
    672657            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
    673  
     658            ! 
    674659            IF(lwp) THEN 
    675660               WRITE(numout,*)  
     
    677662               WRITE(numout,*) '~~~~~~~~~~~~' 
    678663            ENDIF 
    679  
     664            ! 
    680665            ! Update the tracer tendencies 
    681666            DO jk = 1, jpkm1 
     
    700685               ENDIF 
    701686            END DO 
    702  
    703          ENDIF 
    704  
     687            ! 
     688         ENDIF 
     689         ! 
    705690         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work 
    706691            DEALLOCATE( t_bkginc ) 
    707692            DEALLOCATE( s_bkginc ) 
    708693         ENDIF 
    709  
    710  
    711       ELSEIF ( ln_asmdin ) THEN 
    712  
    713          !-------------------------------------------------------------------- 
    714          ! Direct Initialization 
    715          !-------------------------------------------------------------------- 
    716              
     694         !                             !-------------------------------------- 
     695      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization 
     696         !                             !-------------------------------------- 
     697         !             
    717698         IF ( kt == nitdin_r ) THEN 
    718  
     699            ! 
    719700            neuler = 0  ! Force Euler forward step 
    720  
     701            ! 
    721702            ! Initialize the now fields with the background + increment 
    722703            IF (ln_temnofreeze) THEN 
     
    745726!!gm 
    746727 
    747  
    748             IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      & 
    749                &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
    750                &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
    751             IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      & 
    752                &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF) 
    753                &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    754                &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
     728            IF( ln_zps .AND. .NOT. lk_c1d ) THEN      ! Partial steps: before horizontal gradient 
     729               IF(ln_isfcav) THEN                        ! ocean cavities: top and bottom cells (ISF) 
     730                  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,     & 
     731                     &                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     732                     &                     grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) 
     733               ELSE                                      ! no ocean cavities: bottom cells 
     734                  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  !  
     735                     &                        rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
     736               ENDIF 
     737            ENDIF 
    755738 
    756739#if defined key_zdfkpp 
     
    758741!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd 
    759742#endif 
    760  
     743            ! 
    761744            DEALLOCATE( t_bkginc ) 
    762745            DEALLOCATE( s_bkginc ) 
     
    767750      ENDIF 
    768751      ! Perhaps the following call should be in step 
    769       IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment 
     752      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment 
    770753      ! 
    771754   END SUBROUTINE tra_asm_inc 
     
    788771      REAL(wp) :: zincwgt  ! IAU weight for current time step 
    789772      !!---------------------------------------------------------------------- 
    790  
    791       IF ( ln_asmiau ) THEN 
    792  
    793          !-------------------------------------------------------------------- 
    794          ! Incremental Analysis Updating 
    795          !-------------------------------------------------------------------- 
    796  
     773      ! 
     774      !                          !-------------------------------------------- 
     775      IF ( ln_asmiau ) THEN      ! Incremental Analysis Updating 
     776         !                       !-------------------------------------------- 
     777         ! 
    797778         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    798  
     779            ! 
    799780            it = kt - nit000 + 1 
    800781            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
    801  
     782            ! 
    802783            IF(lwp) THEN 
    803784               WRITE(numout,*)  
    804                WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', & 
    805                   &  kt,' with IAU weight = ', wgtiau(it) 
     785               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    806786               WRITE(numout,*) '~~~~~~~~~~~~' 
    807787            ENDIF 
    808  
     788            ! 
    809789            ! Update the dynamic tendencies 
    810790            DO jk = 1, jpkm1 
     
    812792               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 
    813793            END DO 
    814             
     794            ! 
    815795            IF ( kt == nitiaufin_r ) THEN 
    816796               DEALLOCATE( u_bkginc ) 
    817797               DEALLOCATE( v_bkginc ) 
    818798            ENDIF 
    819  
    820          ENDIF 
    821  
    822       ELSEIF ( ln_asmdin ) THEN  
    823  
    824          !-------------------------------------------------------------------- 
    825          ! Direct Initialization 
    826          !-------------------------------------------------------------------- 
    827           
     799            ! 
     800         ENDIF 
     801         !                          !----------------------------------------- 
     802      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization 
     803         !                          !----------------------------------------- 
     804         !          
    828805         IF ( kt == nitdin_r ) THEN 
    829  
     806            ! 
    830807            neuler = 0                    ! Force Euler forward step 
    831  
     808            ! 
    832809            ! Initialize the now fields with the background + increment 
    833810            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 
    834811            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)   
    835  
     812            ! 
    836813            ub(:,:,:) = un(:,:,:)         ! Update before fields 
    837814            vb(:,:,:) = vn(:,:,:) 
    838   
     815            ! 
    839816            DEALLOCATE( u_bkg    ) 
    840817            DEALLOCATE( v_bkg    ) 
     
    864841      REAL(wp) :: zincwgt  ! IAU weight for current time step 
    865842      !!---------------------------------------------------------------------- 
    866  
    867       IF ( ln_asmiau ) THEN 
    868  
    869          !-------------------------------------------------------------------- 
    870          ! Incremental Analysis Updating 
    871          !-------------------------------------------------------------------- 
    872  
     843      ! 
     844      !                             !----------------------------------------- 
     845      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating 
     846         !                          !----------------------------------------- 
     847         ! 
    873848         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    874  
     849            ! 
    875850            it = kt - nit000 + 1 
    876851            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
    877  
     852            ! 
    878853            IF(lwp) THEN 
    879854               WRITE(numout,*)  
     
    882857               WRITE(numout,*) '~~~~~~~~~~~~' 
    883858            ENDIF 
    884  
     859            ! 
    885860            ! Save the tendency associated with the IAU weighted SSH increment 
    886861            ! (applied in dynspg.*) 
     
    891866               DEALLOCATE( ssh_bkginc ) 
    892867            ENDIF 
    893  
    894          ENDIF 
    895  
    896       ELSEIF ( ln_asmdin ) THEN 
    897  
    898          !-------------------------------------------------------------------- 
    899          ! Direct Initialization 
    900          !-------------------------------------------------------------------- 
    901  
     868            ! 
     869         ENDIF 
     870         !                          !----------------------------------------- 
     871      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization 
     872         !                          !----------------------------------------- 
     873         ! 
    902874         IF ( kt == nitdin_r ) THEN 
    903  
    904             neuler = 0                    ! Force Euler forward step 
    905  
    906             ! Initialize the now fields the background + increment 
    907             sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   
    908  
    909             ! Update before fields 
    910             sshb(:,:) = sshn(:,:)          
    911  
     875            ! 
     876            neuler = 0                                   ! Force Euler forward step 
     877            ! 
     878            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment 
     879            ! 
     880            sshb(:,:) = sshn(:,:)                        ! Update before fields 
     881            ! 
    912882            IF( lk_vvl ) THEN 
    913883               DO jk = 1, jpk 
     
    915885               END DO 
    916886            ENDIF 
    917  
     887            ! 
    918888            DEALLOCATE( ssh_bkg    ) 
    919889            DEALLOCATE( ssh_bkginc ) 
    920  
     890            ! 
    921891         ENDIF 
    922892         ! 
     
    937907      !! 
    938908      !!---------------------------------------------------------------------- 
    939       IMPLICIT NONE 
    940       ! 
    941       INTEGER, INTENT(in)           ::   kt   ! Current time step 
     909      INTEGER, INTENT(in)           ::   kt       ! Current time step 
    942910      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
    943911      ! 
     
    949917#endif 
    950918      !!---------------------------------------------------------------------- 
    951  
    952       IF ( ln_asmiau ) THEN 
    953  
    954          !-------------------------------------------------------------------- 
    955          ! Incremental Analysis Updating 
    956          !-------------------------------------------------------------------- 
    957  
     919      ! 
     920      !                             !----------------------------------------- 
     921      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating 
     922         !                          !----------------------------------------- 
     923         ! 
    958924         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
    959  
     925            ! 
    960926            it = kt - nit000 + 1 
    961927            zincwgt = wgtiau(it)      ! IAU weight for the current time step  
    962928            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
    963  
     929            ! 
    964930            IF(lwp) THEN 
    965931               WRITE(numout,*)  
    966                WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 
    967                   &  kt,' with IAU weight = ', wgtiau(it) 
     932               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 
    968933               WRITE(numout,*) '~~~~~~~~~~~~' 
    969934            ENDIF 
    970  
     935            ! 
    971936            ! Sea-ice : LIM-3 case (to add) 
    972  
     937            ! 
    973938#if defined key_lim2 
    974939            ! Sea-ice : LIM-2 case 
     
    1008973 
    1009974#if defined key_cice && defined key_asminc 
    1010             ! Sea-ice : CICE case. Zero ice increment tendency into CICE 
    1011             ndaice_da(:,:) = 0.0_wp 
    1012 #endif 
    1013  
    1014          ENDIF 
    1015  
    1016       ELSEIF ( ln_asmdin ) THEN 
    1017  
    1018          !-------------------------------------------------------------------- 
    1019          ! Direct Initialization 
    1020          !-------------------------------------------------------------------- 
    1021  
     975            ndaice_da(:,:) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE 
     976#endif 
     977 
     978         ENDIF 
     979         !                          !----------------------------------------- 
     980      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization 
     981         !                          !----------------------------------------- 
     982         ! 
    1022983         IF ( kt == nitdin_r ) THEN 
    1023  
     984            ! 
    1024985            neuler = 0                    ! Force Euler forward step 
    1025  
     986            ! 
    1026987            ! Sea-ice : LIM-3 case (to add) 
    1027  
     988            ! 
    1028989#if defined key_lim2 
    1029990            ! Sea-ice : LIM-2 case. 
     
    10411002               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
    10421003            ELSEWHERE 
    1043                zhicifinc(:,:) = 0.0_wp 
     1004               zhicifinc(:,:) = 0._wp 
    10441005            END WHERE 
    10451006            ! 
     
    10501011            ! seaice salinity balancing (to add) 
    10511012#endif 
    1052   
     1013            ! 
    10531014#if defined key_cice && defined key_asminc 
    10541015            ! Sea-ice : CICE case. Pass ice increment tendency into CICE 
    10551016           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 
    10561017#endif 
    1057            IF ( .NOT. PRESENT(kindic) ) THEN 
    1058               DEALLOCATE( seaice_bkginc ) 
    1059            END IF 
    1060  
     1018            IF ( .NOT. PRESENT(kindic) ) THEN 
     1019               DEALLOCATE( seaice_bkginc ) 
     1020            END IF 
     1021            ! 
    10611022         ELSE 
    1062  
     1023            ! 
    10631024#if defined key_cice && defined key_asminc 
    1064             ! Sea-ice : CICE case. Zero ice increment tendency into CICE  
    1065             ndaice_da(:,:) = 0.0_wp 
    1066 #endif 
    1067           
     1025            ndaice_da(:,:) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE 
     1026 
     1027#endif 
     1028            ! 
    10681029         ENDIF 
    10691030 
     
    11421103! 
    11431104!#endif 
    1144  
     1105         ! 
    11451106      ENDIF 
    1146  
     1107      ! 
    11471108   END SUBROUTINE seaice_asm_inc 
    11481109    
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90

    r2287 r5777  
    66 
    77   IMPLICIT NONE 
    8  
    9    !! * Routine accessibility 
    108   PRIVATE 
    119 
    12    !! * Shared Modules variables 
    13    CHARACTER (LEN=40), PUBLIC, PARAMETER :: & 
    14       & c_asmbkg = 'assim_background_state_Jb',  & !: Filename for storing the  
    15                                                    !: background state for use  
    16                                                    !: in the Jb term 
    17       & c_asmdin = 'assim_background_state_DI',  & !: Filename for storing the  
    18                                                    !: background state for direct  
    19                                                    !: initialization 
    20       & c_asmtrj = 'assim_trj',                  & !: Filename for storing the  
    21                                                    !: reference trajectory 
    22       & c_asminc = 'assim_background_increments'   !: Filename for storing the  
    23                                                    !: increments to the background 
    24                                                    !: state 
     10   CHARACTER(LEN=40), PUBLIC, PARAMETER ::   c_asmbkg = 'assim_background_state_Jb'   !: Filename for storing the background state 
     11   !                                                                                  !  for use in the Jb term 
     12   CHARACTER(LEN=40), PUBLIC, PARAMETER ::   c_asmdin = 'assim_background_state_DI'   !: Filename for storing the background state 
     13   !                                                                                  !  for direct initialization 
     14   CHARACTER(LEN=40), PUBLIC, PARAMETER ::   c_asmtrj = 'assim_trj'                   !: Filename for storing the reference trajectory 
     15   CHARACTER(LEN=40), PUBLIC, PARAMETER ::   c_asminc = 'assim_background_increments' !: Filename for storing the increments  
     16   !                                                                                  !  to the background state 
    2517 
    26    INTEGER, PUBLIC :: nitbkg_r      !: Background time step referenced to nit000 
    27    INTEGER, PUBLIC :: nitdin_r      !: Direct Initialization time step referenced to nit000 
    28    INTEGER, PUBLIC :: nitiaustr_r   !: IAU starting time step referenced to nit000 
    29    INTEGER, PUBLIC :: nitiaufin_r   !: IAU final time step referenced to nit000 
    30    INTEGER, PUBLIC :: nittrjfrq     !: Frequency of trajectory output for 4D-VAR 
     18   INTEGER, PUBLIC ::   nitbkg_r      !: Background time step referenced to nit000 
     19   INTEGER, PUBLIC ::   nitdin_r      !: Direct Initialization time step referenced to nit000 
     20   INTEGER, PUBLIC ::   nitiaustr_r   !: IAU starting time step referenced to nit000 
     21   INTEGER, PUBLIC ::   nitiaufin_r   !: IAU final time step referenced to nit000 
     22   INTEGER, PUBLIC ::   nittrjfrq     !: Frequency of trajectory output for 4D-VAR 
    3123 
    3224   !!---------------------------------------------------------------------- 
     
    3426   !! $Id$ 
    3527   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    36    !!---------------------------------------------------------------------- 
    37  
     28   !!====================================================================== 
    3829END MODULE asmpar 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r4699 r5777  
    88   !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions 
    99   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    10    !!            3.6  !  2012-01  (C. Rousset) add ice boundary conditions for lim3 
     10   !!            3.6  !  2014-01  (C. Rousset) add ice boundary conditions for lim3 
    1111   !!---------------------------------------------------------------------- 
    1212#if defined key_bdy  
     
    2222 
    2323   TYPE, PUBLIC ::   OBC_INDEX    !: Indices and weights which define the open boundary 
    24       INTEGER,          DIMENSION(jpbgrd) ::  nblen 
    25       INTEGER,          DIMENSION(jpbgrd) ::  nblenrim 
    26       INTEGER, POINTER, DIMENSION(:,:)   ::  nbi 
    27       INTEGER, POINTER, DIMENSION(:,:)   ::  nbj 
    28       INTEGER, POINTER, DIMENSION(:,:)   ::  nbr 
    29       INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
    30       REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbw 
    31       REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbd 
    32       REAL(wp)   , POINTER, DIMENSION(:,:)   ::  nbdout 
    33       REAL(wp)   , POINTER, DIMENSION(:,:)   ::  flagu 
    34       REAL(wp)   , POINTER, DIMENSION(:,:)   ::  flagv 
     24      INTEGER ,          DIMENSION(jpbgrd) ::  nblen 
     25      INTEGER ,          DIMENSION(jpbgrd) ::  nblenrim 
     26      INTEGER , POINTER, DIMENSION(:,:)    ::  nbi 
     27      INTEGER , POINTER, DIMENSION(:,:)    ::  nbj 
     28      INTEGER , POINTER, DIMENSION(:,:)    ::  nbr 
     29      INTEGER , POINTER, DIMENSION(:,:)    ::  nbmap 
     30      REAL(wp), POINTER, DIMENSION(:,:)    ::  nbw 
     31      REAL(wp), POINTER, DIMENSION(:,:)    ::  nbd 
     32      REAL(wp), POINTER, DIMENSION(:,:)    ::  nbdout 
     33      REAL(wp), POINTER, DIMENSION(:,:)    ::  flagu 
     34      REAL(wp), POINTER, DIMENSION(:,:)    ::  flagv 
    3535   END TYPE OBC_INDEX 
    3636 
     
    4141 
    4242   TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data 
    43       INTEGER,       DIMENSION(2)     ::  nread 
    44       LOGICAL                         ::  ll_ssh 
    45       LOGICAL                         ::  ll_u2d 
    46       LOGICAL                         ::  ll_v2d 
    47       LOGICAL                         ::  ll_u3d 
    48       LOGICAL                         ::  ll_v3d 
    49       LOGICAL                         ::  ll_tem 
    50       LOGICAL                         ::  ll_sal 
    51       REAL(wp), POINTER, DIMENSION(:)     ::  ssh 
    52       REAL(wp), POINTER, DIMENSION(:)     ::  u2d 
    53       REAL(wp), POINTER, DIMENSION(:)     ::  v2d 
    54       REAL(wp), POINTER, DIMENSION(:,:)   ::  u3d 
    55       REAL(wp), POINTER, DIMENSION(:,:)   ::  v3d 
    56       REAL(wp), POINTER, DIMENSION(:,:)   ::  tem 
    57       REAL(wp), POINTER, DIMENSION(:,:)   ::  sal 
     43      INTEGER          , DIMENSION(2)   ::  nread 
     44      LOGICAL                           ::  ll_ssh 
     45      LOGICAL                           ::  ll_u2d 
     46      LOGICAL                           ::  ll_v2d 
     47      LOGICAL                           ::  ll_u3d 
     48      LOGICAL                           ::  ll_v3d 
     49      LOGICAL                           ::  ll_tem 
     50      LOGICAL                           ::  ll_sal 
     51      REAL(wp), POINTER, DIMENSION(:)   ::  ssh 
     52      REAL(wp), POINTER, DIMENSION(:)   ::  u2d 
     53      REAL(wp), POINTER, DIMENSION(:)   ::  v2d 
     54      REAL(wp), POINTER, DIMENSION(:,:) ::  u3d 
     55      REAL(wp), POINTER, DIMENSION(:,:) ::  v3d 
     56      REAL(wp), POINTER, DIMENSION(:,:) ::  tem 
     57      REAL(wp), POINTER, DIMENSION(:,:) ::  sal 
    5858#if defined key_lim2 
    59       LOGICAL                         ::  ll_frld 
    60       LOGICAL                         ::  ll_hicif 
    61       LOGICAL                         ::  ll_hsnif 
    62       REAL(wp), POINTER, DIMENSION(:)     ::  frld 
    63       REAL(wp), POINTER, DIMENSION(:)     ::  hicif 
    64       REAL(wp), POINTER, DIMENSION(:)     ::  hsnif 
     59      LOGICAL                           ::   ll_frld 
     60      LOGICAL                           ::   ll_hicif 
     61      LOGICAL                           ::   ll_hsnif 
     62      REAL(wp), POINTER, DIMENSION(:)   ::   frld 
     63      REAL(wp), POINTER, DIMENSION(:)   ::   hicif 
     64      REAL(wp), POINTER, DIMENSION(:)   ::   hsnif 
    6565#elif defined key_lim3 
    66       LOGICAL                         ::  ll_a_i 
    67       LOGICAL                         ::  ll_ht_i 
    68       LOGICAL                         ::  ll_ht_s 
    69       REAL, POINTER, DIMENSION(:,:)   ::  a_i   !: now ice leads fraction climatology 
    70       REAL, POINTER, DIMENSION(:,:)   ::  ht_i  !: Now ice  thickness climatology 
    71       REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness 
     66      LOGICAL                           ::   ll_a_i 
     67      LOGICAL                           ::   ll_ht_i 
     68      LOGICAL                           ::   ll_ht_s 
     69      REAL(wp), POINTER, DIMENSION(:,:) ::   a_i    !: now ice leads fraction climatology 
     70      REAL(wp), POINTER, DIMENSION(:,:) ::   ht_i   !: Now ice  thickness climatology 
     71      REAL(wp), POINTER, DIMENSION(:,:) ::   ht_s   !: now snow thickness 
    7272#endif 
    7373   END TYPE OBC_DATA 
     
    9999   INTEGER, DIMENSION(jp_bdy)           ::   nn_tra_dta     !: = 0 use the initial state as bdy dta ;  
    100100                                                            !: = 1 read it in a NetCDF file 
    101    LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    102    LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
    103    REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    104    REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 
     101   LOGICAL , DIMENSION(jp_bdy) ::   ln_tra_dmp              !: =T Tracer damping 
     102   LOGICAL , DIMENSION(jp_bdy) ::   ln_dyn3d_dmp            !: =T Baroclinic velocity damping 
     103   REAL(wp), DIMENSION(jp_bdy) ::   rn_time_dmp             !: Damping time scale in days 
     104   REAL(wp), DIMENSION(jp_bdy) ::   rn_time_dmp_out         !: Damping time scale in days at radiation outflow points 
    105105 
    106106   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_ice_lim       ! Choice of boundary condition for sea ice variables  
    107    INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim_dta   !: = 0 use the initial state as bdy dta ;  
     107   INTEGER , DIMENSION(jp_bdy)          ::   nn_ice_lim_dta   !: = 0 use the initial state as bdy dta ;  
    108108                                                              !: = 1 read it in a NetCDF file 
    109    REAL(wp),    DIMENSION(jp_bdy) ::   rn_ice_tem             !: choice of the temperature of incoming sea ice 
    110    REAL(wp),    DIMENSION(jp_bdy) ::   rn_ice_sal             !: choice of the salinity    of incoming sea ice 
    111    REAL(wp),    DIMENSION(jp_bdy) ::   rn_ice_age             !: choice of the age         of incoming sea ice 
     109   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_tem              !: choice of the temperature of incoming sea ice 
     110   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_sal              !: choice of the salinity    of incoming sea ice 
     111   REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_age              !: choice of the age         of incoming sea ice 
    112112   ! 
    113113    
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r5656 r5777  
    5959      !! 
    6060      !!---------------------------------------------------------------------- 
    61       INTEGER, INTENT( in ) :: kt     ! Main time step counter 
    62       INTEGER               :: ib_bdy ! Loop index 
    63  
     61      INTEGER, INTENT( in ) ::   kt   ! Main time step counter 
     62      ! 
     63      INTEGER ::   ib_bdy   ! Loop index 
     64      !!---------------------------------------------------------------------- 
     65      ! 
    6466#if defined key_lim3 
    6567      CALL lim_var_glo2eqv 
    6668#endif 
    67  
     69      ! 
    6870      DO ib_bdy=1, nb_bdy 
    69  
     71         ! 
    7072         SELECT CASE( cn_ice_lim(ib_bdy) ) 
    7173         CASE('none') 
     
    7678            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' ) 
    7779         END SELECT 
    78  
     80         ! 
    7981      END DO 
    80  
     82      ! 
    8183#if defined key_lim3 
    8284      CALL lim_var_zapsmall 
    8385      CALL lim_var_agg(1) 
    8486#endif 
    85  
     87      ! 
    8688   END SUBROUTINE bdy_ice_lim 
     89 
    8790 
    8891   SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy ) 
     
    9699      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382. 
    97100      !!------------------------------------------------------------------------------ 
    98       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    99       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    100       INTEGER,         INTENT(in) ::   kt   ! main time-step counter 
     101      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices 
     102      TYPE(OBC_DATA),  INTENT(in) ::   dta     ! OBC external data 
     103      INTEGER,         INTENT(in) ::   kt      ! main time-step counter 
    101104      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    102  
     105      ! 
    103106      INTEGER  ::   jpbound            ! 0 = incoming ice 
    104                                        ! 1 = outgoing ice 
     107      !                                ! 1 = outgoing ice 
    105108      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices 
    106109      INTEGER  ::   ji, jj, ii, ij     ! local scalar 
     
    111114     USE ice_2, vt_i => hicm 
    112115#endif 
    113  
    114       !!------------------------------------------------------------------------------ 
    115       ! 
    116       IF( nn_timing == 1 ) CALL timing_start('bdy_ice_frs') 
     116      !!------------------------------------------------------------------------------ 
     117      ! 
     118      IF( nn_timing == 1 )   CALL timing_start('bdy_ice_frs') 
    117119      ! 
    118120      jgrd = 1      ! Everything is at T-points here 
     
    181183            ! condition on ice thickness depends on the ice velocity 
    182184            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values 
    183             jpbound = 0; ii = ji; ij = jj; 
    184  
     185            jpbound = 0   ;   ii = ji   ;   ij = jj 
     186            ! 
    185187            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj 
    186188            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj 
    187189            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1 
    188190            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1 
    189  
     191            ! 
    190192            IF( nn_ice_lim_dta(ib_bdy) == 0 ) jpbound = 0; ii = ji; ij = jj   ! case ice boundaries = initial conditions 
    191                                                                               !      do not make state variables dependent on velocity 
    192                 
    193  
     193            !                                                                 !      do not make state variables dependent on velocity 
     194            ! 
    194195            rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice 
    195  
     196            ! 
    196197            ! concentration and thickness 
    197198            a_i (ji,jj,jl) = a_i (ii,ij,jl) * rswitch 
    198199            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * rswitch 
    199200            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * rswitch 
    200  
     201            ! 
    201202            ! Ice and snow volumes 
    202203            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    203204            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl) 
    204  
     205            ! 
    205206            SELECT CASE( jpbound ) 
    206  
    207             CASE( 0 ) ! velocity is inward 
    208  
     207            ! 
     208            CASE( 0 )   ! velocity is inward 
     209               ! 
    209210               ! Ice salinity, age, temperature 
    210211               sm_i(ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * rn_simin 
     
    218219                  s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin 
    219220               END DO 
    220                 
    221             CASE( 1 ) ! velocity is outward 
    222   
     221               ! 
     222            CASE( 1 )   ! velocity is outward 
     223               ! 
    223224               ! Ice salinity, age, temperature 
    224225               sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * rn_simin 
     
    232233                  s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin 
    233234               END DO 
    234  
     235               ! 
    235236            END SELECT 
    236  
    237             ! if salinity is constant, then overwrite rn_ice_sal 
    238             IF( nn_icesal == 1 ) THEN 
    239                sm_i(ji,jj,jl)   = rn_icesal 
     237            ! 
     238            IF( nn_icesal == 1 ) THEN     ! constant salinity : overwrite rn_ice_sal 
     239               sm_i(ji,jj  ,jl) = rn_icesal 
    240240               s_i (ji,jj,:,jl) = rn_icesal 
    241241            ENDIF 
    242  
     242            ! 
    243243            ! contents 
    244244            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
     
    259259               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) * r1_nlay_i 
    260260            END DO 
    261  
     261            ! 
    262262         END DO 
    263   
     263         ! 
    264264         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy ) 
    265265         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
     
    267267         CALL lbc_bdy_lnk(  v_i(:,:,jl), 'T', 1., ib_bdy ) 
    268268         CALL lbc_bdy_lnk(  v_s(:,:,jl), 'T', 1., ib_bdy ) 
    269   
     269         ! 
    270270         CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy ) 
    271271         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy ) 
     
    280280            CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy ) 
    281281         END DO 
    282  
     282         ! 
    283283      END DO !jl 
    284      
     284      ! 
    285285#endif 
    286286      !       
    287       IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs') 
     287      IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs') 
    288288      ! 
    289289   END SUBROUTINE bdy_ice_frs 
     
    300300      !! 2013-06 : C. Rousset 
    301301      !!------------------------------------------------------------------------------ 
    302       !! 
    303302      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
     303      ! 
    304304      INTEGER  ::   jb, jgrd           ! dummy loop indices 
    305305      INTEGER  ::   ji, jj             ! local scalar 
    306306      INTEGER  ::   ib_bdy             ! Loop index 
    307307      REAL(wp) ::   zmsk1, zmsk2, zflag 
    308      !!------------------------------------------------------------------------------ 
     308      !!------------------------------------------------------------------------------ 
    309309      ! 
    310310      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_lim_dyn') 
     
    313313         ! 
    314314         SELECT CASE( cn_ice_lim(ib_bdy) ) 
    315  
     315         ! 
    316316         CASE('none') 
    317  
    318317            CYCLE 
    319              
     318            ! 
    320319         CASE('frs') 
    321              
     320            ! 
    322321            IF( nn_ice_lim_dta(ib_bdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions  
    323                                                                !      do not change ice velocity (it is only computed by rheology) 
    324   
     322            !                                                  !      do not change ice velocity (it is only computed by rheology) 
    325323            SELECT CASE ( cd_type ) 
    326                 
    327             CASE ( 'U' ) 
    328                 
     324            !      
     325            CASE ( 'U' )   
    329326               jgrd = 2      ! u velocity 
    330327               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
     
    332329                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    333330                  zflag = idx_bdy(ib_bdy)%flagu(jb,jgrd) 
    334                    
     331                  ! 
    335332                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
    336333                     ! one of the two zmsk is always 0 (because of zflag) 
    337334                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice 
    338335                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice 
    339                       
     336                      
    340337                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    341338                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
     
    349346                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice 
    350347                  u_ice(ji,jj) = rswitch * u_ice(ji,jj) 
    351                    
    352                ENDDO 
    353                 
     348                  ! 
     349               END DO 
    354350               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) 
    355                 
     351               ! 
    356352            CASE ( 'V' ) 
    357                 
    358353               jgrd = 3      ! v velocity 
    359354               DO jb = 1, idx_bdy(ib_bdy)%nblen(jgrd) 
     
    361356                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    362357                  zflag = idx_bdy(ib_bdy)%flagv(jb,jgrd) 
    363                    
     358                  ! 
    364359                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries 
    365360                     ! one of the two zmsk is always 0 (because of zflag) 
    366361                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice 
    367362                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice 
    368                       
     363                      
    369364                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    370365                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
     
    378373                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice 
    379374                  v_ice(ji,jj) = rswitch * v_ice(ji,jj) 
    380                    
    381                ENDDO 
    382                 
     375                  ! 
     376               END DO 
    383377               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) 
    384                    
     378               ! 
    385379            END SELECT 
    386              
     380            ! 
    387381         CASE DEFAULT 
    388382            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' ) 
    389383         END SELECT 
    390           
    391       ENDDO 
    392  
     384         ! 
     385      END DO 
     386      ! 
    393387      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn') 
    394        
     388      ! 
    395389    END SUBROUTINE bdy_ice_lim_dyn 
    396390 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r5656 r5777  
    7676      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 
    7777      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers 
    78       INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy         !   -       - 
     78      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       - 
    7979      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
    8080      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r5737 r5777  
    1515   !!   'key_dynspg_flt'                              filtered free surface 
    1616   !!---------------------------------------------------------------------- 
    17    USE timing          ! Timing 
    1817   USE oce             ! ocean dynamics and tracers  
    19    USE sbcisf          ! ice shelf 
     18   USE bdy_oce         ! ocean open boundary conditions 
     19   USE sbc_oce         ! ocean surface boundary conditions 
    2020   USE dom_oce         ! ocean space and time domain  
    2121   USE phycst          ! physical constants 
    22    USE bdy_oce         ! ocean open boundary conditions 
     22   USE sbcisf          ! ice shelf 
     23   ! 
     24   USE in_out_manager  ! I/O manager 
    2325   USE lib_mpp         ! for mppsum 
    24    USE in_out_manager  ! I/O manager 
    25    USE sbc_oce         ! ocean surface boundary conditions 
     26   USE timing          ! Timing 
     27   USE lib_fortran     ! Fortran routines library 
    2628 
    2729   IMPLICIT NONE 
     
    3335#  include "domzgr_substitute.h90" 
    3436   !!---------------------------------------------------------------------- 
    35    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     37   !! NEMO/OPA 3.6 , NEMO Consortium (2014) 
    3638   !! $Id$  
    3739   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     
    7880      TYPE(OBC_INDEX), POINTER :: idx 
    7981      !!----------------------------------------------------------------------------- 
    80  
    81       IF( nn_timing == 1 ) CALL timing_start('bdy_vol') 
    82  
     82      ! 
     83      IF( nn_timing == 1 )   CALL timing_start('bdy_vol') 
     84      ! 
    8385      IF( ln_vol ) THEN 
    84  
     86      ! 
    8587      IF( kt == nit000 ) THEN  
    8688         IF(lwp) WRITE(numout,*) 
     
    9193      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
    9294      ! ----------------------------------------------------------------------- 
     95!!gm replace these lines : 
    9396      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 
    9497      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain 
     98!!gm   by : 
     99!!gm      z_cflxemp = glob_sum(  ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
     100!!gm 
    95101 
    96102      ! Transport through the unstructured open boundary 
    97103      ! ------------------------------------------------ 
    98       zubtpecor = 0.e0 
     104      zubtpecor = 0._wp 
    99105      DO ib_bdy = 1, nb_bdy 
    100106         idx => idx_bdy(ib_bdy) 
    101  
     107         ! 
    102108         jgrd = 2                               ! cumulate u component contribution first  
    103109         DO jb = 1, idx%nblenrim(jgrd) 
     
    116122            END DO 
    117123         END DO 
    118  
     124         ! 
    119125      END DO 
    120126      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain 
     
    123129      ! ------------------------------ 
    124130      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot  
    125       ELSE                   ;   zubtpecor =   zubtpecor             / bdysurftot 
     131      ELSE                      ;   zubtpecor =   zubtpecor             / bdysurftot 
    126132      END IF 
    127133 
    128134      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
    129135      ! ------------------------------------------------------------- 
    130       ztranst = 0.e0 
     136      ztranst = 0._wp 
    131137      DO ib_bdy = 1, nb_bdy 
    132138         idx => idx_bdy(ib_bdy) 
    133  
     139         ! 
    134140         jgrd = 2                               ! correct u component 
    135141         DO jb = 1, idx%nblenrim(jgrd) 
     
    150156            END DO 
    151157         END DO 
    152  
     158         ! 
    153159      END DO 
    154160      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain 
     
    169175      ! 
    170176      END IF ! ln_vol 
    171  
     177      ! 
    172178   END SUBROUTINE bdy_vol 
    173179 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/C1D/c1d.F90

    r5215 r5777  
    4848      !!---------------------------------------------------------------------- 
    4949      ! 
    50  
    5150      REWIND( numnam_ref )              ! Namelist namc1d in reference namelist : Tracer advection scheme 
    5251      READ  ( numnam_ref, namc1d, IOSTAT = ios, ERR = 901) 
     
    5756902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namc1d in configuration namelist', lwp ) 
    5857      IF(lwm) WRITE ( numond, namc1d ) 
    59  
    6058      ! 
    6159      IF(lwp) THEN                    ! Control print 
     
    6967      ENDIF 
    7068      ! 
    71       ! 
    7269   END SUBROUTINE c1d_init 
    7370 
     
    7774   !!---------------------------------------------------------------------- 
    7875   USE par_kind         ! kind parameters 
    79  
    8076   LOGICAL, PUBLIC, PARAMETER ::   lk_c1d = .FALSE.   !: 1D config. flag de-activated 
    8177   REAL(wp)                   ::   rn_lat1d, rn_lon1d 
    8278   LOGICAL , PUBLIC           ::   ln_c1d_locpt = .FALSE.  
    83  
    8479CONTAINS 
    85  
    8680   SUBROUTINE c1d_init               ! Dummy routine 
    8781   END SUBROUTINE c1d_init 
    88  
    8982#endif 
    9083 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90

    r5758 r5777  
    165165      CALL iom_put( "eken", rke_crs ) 
    166166 
    167       !  Horizontal divergence ( following OPA_SRC/DYN/divcur.F90 )  
     167      !  Horizontal divergence ( following OPA_SRC/DYN/divhor.F90 )  
    168168      DO jk = 1, jpkm1 
    169169         DO ji = 2, jpi_crsm1 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5737 r5777  
    77   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 
    88   !!---------------------------------------------------------------------- 
    9 #if defined key_diaar5   || defined key_esopa 
     9#if defined key_diaar5 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_diaar5'  :                           activate ar5 diagnotics 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r5363 r5777  
    88   !!                 !  1997-08  (G. Madec)  optimization 
    99   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content  
    10    !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module 
    11    !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
    12    !!---------------------------------------------------------------------- 
    13 #if   defined key_diahth   ||   defined key_esopa 
     10   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     11   !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
     12   !!---------------------------------------------------------------------- 
     13#if   defined key_diahth 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   'key_diahth' :                              thermocline depth diag. 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r5758 r5777  
    2929   USE dynadv, ONLY: ln_dynadv_vec 
    3030   USE zdf_oce         ! ocean vertical physics 
    31    USE ldftra          ! ocean active tracers: lateral physics 
    32    USE ldfdyn_oce      ! ocean dynamics: lateral physics 
     31   USE ldftra          ! lateral physics: eddy diffusivity coef. 
    3332   USE sol_oce         ! solver variables 
    3433   USE sbc_oce         ! Surface boundary condition: ocean fields 
     
    402401      !!      Each nwrite time step, output the instantaneous or mean fields 
    403402      !!---------------------------------------------------------------------- 
    404       !! 
    405       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    406       !! 
     403      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     404      ! 
    407405      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout 
    408406      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names 
     
    872870      !!---------------------------------------------------------------------- 
    873871      !  
    874 !     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep 
    875  
    876872      ! 0. Initialisation 
    877873      ! ----------------- 
     
    974970      ENDIF 
    975971#endif 
    976         
    977 !     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep 
    978972      !  
    979973   END SUBROUTINE dia_wri_state 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r5506 r5777  
    7171      !!                                   =2 put at location runoff 
    7272      !!---------------------------------------------------------------------- 
    73       INTEGER ::   jc            ! dummy loop indices 
    74       INTEGER :: isrow           ! local index 
    75       !!---------------------------------------------------------------------- 
    76        
     73      INTEGER ::   jc      ! dummy loop indices 
     74      INTEGER ::   isrow   ! local index 
     75      !!---------------------------------------------------------------------- 
     76      ! 
    7777      IF(lwp) WRITE(numout,*) 
    7878      IF(lwp) WRITE(numout,*)'dom_clo : closed seas ' 
    7979      IF(lwp) WRITE(numout,*)'~~~~~~~' 
    80  
     80      ! 
    8181      ! initial values 
    8282      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1 
    8383      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1 
    84  
     84      ! 
    8585      ! set the closed seas (in data domain indices) 
    8686      ! ------------------- 
    87  
     87      ! 
    8888      IF( cp_cfg == "orca" ) THEN 
    8989         ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5770 r5777  
    1919   !!   dom_nam        : read and contral domain namelists 
    2020   !!   dom_ctl        : control print for the ocean domain 
     21   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    2122   !!---------------------------------------------------------------------- 
    2223   USE oce             ! ocean variables 
     
    2526   USE phycst          ! physical constants 
    2627   USE closea          ! closed seas 
    27    USE in_out_manager  ! I/O manager 
    28    USE lib_mpp         ! distributed memory computing library 
    29  
    3028   USE domhgr          ! domain: set the horizontal mesh 
    3129   USE domzgr          ! domain: set the vertical mesh 
     
    3634   USE c1d             ! 1D vertical configuration 
    3735   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
     36   ! 
     37   USE in_out_manager  ! I/O manager 
     38   USE lib_mpp         ! distributed memory computing library 
     39   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    3840   USE timing          ! Timing 
    39    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    4041 
    4142   IMPLICIT NONE 
     
    8889      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
    8990      ! 
    90       ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points 
    91       hu_0(:,:) = 0.0_wp                       ! Reference ocean depth at U-points 
    92       hv_0(:,:) = 0.0_wp                       ! Reference ocean depth at V-points 
     91      ht_0(:,:) = 0._wp                        ! Reference ocean depth at T-points 
     92      hu_0(:,:) = 0._wp                        ! Reference ocean depth at U-points 
     93      hv_0(:,:) = 0._wp                        ! Reference ocean depth at V-points 
    9394      DO jk = 1, jpk 
    9495         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
     
    9798      END DO 
    9899      ! 
    99       IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
     100      IF( lk_vvl         )   CALL dom_vvl_init ! Vertical variable mesh 
    100101      ! 
    101102      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 
     
    153154      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run 
    154155      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901) 
    155 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in reference namelist', lwp ) 
     156901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp ) 
    156157 
    157158      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run 
    158159      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 
    159 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
     160902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
    160161      IF(lwm) WRITE ( numond, namrun ) 
    161162      ! 
     
    249250904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    250251      IF(lwm) WRITE ( numond, namdom ) 
    251  
     252      ! 
    252253      IF(lwp) THEN 
    253254         WRITE(numout,*) 
     
    291292         WRITE(numout,*) '                                      ppacr2            = ', ppacr2 
    292293      ENDIF 
    293  
     294      ! 
    294295      ntopo     = nn_bathy          ! conversion DOCTOR names into model names (this should disappear soon) 
    295296      e3zps_min = rn_e3zps_min 
     
    385386   END SUBROUTINE dom_ctl 
    386387 
     388 
    387389   SUBROUTINE dom_stiff 
    388390      !!---------------------------------------------------------------------- 
     
    403405      REAL(wp), DIMENSION(4) :: zr1 
    404406      !!---------------------------------------------------------------------- 
    405       rx1(:,:) = 0.e0 
    406       zrxmax   = 0.e0 
    407       zr1(:)   = 0.e0 
    408        
     407      rx1(:,:) = 0._wp 
     408      zrxmax   = 0._wp 
     409      zr1(:)   = 0._wp 
     410      ! 
    409411      DO ji = 2, jpim1 
    410412         DO jj = 2, jpjm1 
     
    431433         END DO 
    432434      END DO 
    433  
    434435      CALL lbc_lnk( rx1, 'T', 1. ) 
    435  
    436       zrxmax = MAXVAL(rx1) 
    437  
     436      ! 
     437      zrxmax = MAXVAL( rx1 ) 
     438      ! 
    438439      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain 
    439  
     440      ! 
    440441      IF(lwp) THEN 
    441442         WRITE(numout,*) 
     
    443444         WRITE(numout,*) '~~~~~~~~~' 
    444445      ENDIF 
    445  
     446      ! 
    446447   END SUBROUTINE dom_stiff 
    447  
    448  
    449448 
    450449   !!====================================================================== 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5770 r5777  
    77   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
    88   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file 
    9    !!            3.0  ! 2008-01  (S. Masson) add dom_uniq  
    10    !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
     9   !!            3.0  ! 2008-01  (S. Masson)  add dom_uniq  
    1110   !!---------------------------------------------------------------------- 
    1211 
    1312   !!---------------------------------------------------------------------- 
    1413   !!   dom_wri        : create and write mesh and mask file(s) 
    15    !!   dom_uniq       : 
     14   !!   dom_uniq       : identify unique point of a grid (TUVF) 
    1615   !!---------------------------------------------------------------------- 
    1716   USE dom_oce         ! ocean space and time domain 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5770 r5777  
    2929   USE daymod          ! calendar 
    3030   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine) 
    31    USE ldftra          ! ocean active tracers: lateral physics 
     31   USE ldftra          ! lateral physics: ocean active tracers 
    3232   USE zdf_oce         ! ocean vertical physics 
    3333   USE phycst          ! physical constants 
     
    7474      ! 
    7575 
    76       IF(lwp) WRITE(numout,*) ' ' 
     76      IF(lwp) WRITE(numout,*) 
    7777      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 
    7878      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    7979 
    80       CALL dta_tsd_init                       ! Initialisation of T & S input data 
    81       IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data 
     80                     CALL dta_tsd_init        ! Initialisation of T & S input data 
     81      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data 
    8282 
    8383      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
     
    101101         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp 
    102102         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp   
    103          rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp 
    104          hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp 
     103                                    hdivn(:,:,:) = 0._wp 
    105104         ! 
    106105         IF( cp_cfg == 'eel' ) THEN 
     
    150149      ! 
    151150      ! 
    152       un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 
    153       ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 
     151      un_b(:,:) = 0._wp   ;  vn_b(:,:) = 0._wp 
     152      ub_b(:,:) = 0._wp   ;  vb_b(:,:) = 0._wp 
    154153      ! 
    155154      DO jk = 1, jpkm1 
     
    189188      !! References :  Philander ??? 
    190189      !!---------------------------------------------------------------------- 
    191       INTEGER  :: ji, jj, jk 
    192       REAL(wp) ::   zsal = 35.50 
     190      INTEGER  ::   ji, jj, jk 
     191      REAL(wp) ::   zsal = 35.50_wp 
    193192      !!---------------------------------------------------------------------- 
    194193      ! 
     
    220219      !!                and relative vorticity fields 
    221220      !!---------------------------------------------------------------------- 
    222       USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine) 
     221      USE divhor     ! hor. divergence      (div_hor routine) 
    223222      USE iom 
    224223      ! 
     
    269268            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
    270269            ! 
    271             ! set the dynamics: U,V, hdiv, rot (and ssh if necessary) 
     270            ! set the dynamics: U,V, hdiv (and ssh if necessary) 
    272271            ! ---------------- 
    273272            ! Start EEL5 configuration with barotropic geostrophic velocities  
     
    305304            ENDIF 
    306305            ! 
    307             CALL div_cur( nit000 )                  ! horizontal divergence and relative vorticity (curl) 
     306!!gm  Check  here call to div_hor should not be necessary 
     307!!gm         div_hor call runoffs  not sure they are defined at that level 
     308            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl) 
    308309            ! N.B. the vertical velocity will be computed from the horizontal divergence field 
    309310            ! in istate by a call to wzv routine 
     
    358359      !! 
    359360      !! ** Method  : - set temprature field 
    360       !!              - set salinity field 
     361      !!              - set salinity   field 
    361362      !!---------------------------------------------------------------------- 
    362363      INTEGER :: ji, jj, jk  ! dummy loop indices 
     
    445446      !!---------------------------------------------------------------------- 
    446447      USE dynspg          ! surface pressure gradient             (dyn_spg routine) 
    447       USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine) 
     448      USE divhor          ! hor. divergence                       (div_hor routine) 
    448449      USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    449450      ! 
     
    531532      indic = 0 
    532533      CALL dyn_spg( nit000, indic )       ! surface pressure gradient 
    533  
     534      ! 
    534535      ! the new velocity is ua*rdt 
    535  
     536      ! 
    536537      CALL lbc_lnk( ua, 'U', -1. ) 
    537538      CALL lbc_lnk( va, 'V', -1. ) 
     
    543544      un(:,:,:) = ub(:,:,:) 
    544545      vn(:,:,:) = vb(:,:,:) 
    545         
    546       ! Compute the divergence and curl 
    547  
    548       CALL div_cur( nit000 )            ! now horizontal divergence and curl 
    549  
    550       hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value 
    551       rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value 
     546      ! 
     547!!gm  Check  here call to div_hor should not be necessary 
     548!!gm         div_hor call runoffs  not sure they are defined at that level 
     549      CALL div_hor( nit000 )            ! now horizontal divergence 
    552550      ! 
    553551      CALL wrk_dealloc( jpi,jpj,jpk,   zprn) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90

    r5775 r5777  
    1 MODULE divcur 
     1MODULE divhor 
    22   !!============================================================================== 
    3    !!                       ***  MODULE  divcur  *** 
    4    !! Ocean diagnostic variable : horizontal divergence and relative vorticity 
     3   !!                       ***  MODULE  divhor  *** 
     4   !! Ocean diagnostic variable : now horizontal divergence 
    55   !!============================================================================== 
    6    !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code 
    7    !!            4.0  ! 1991-11  (G. Madec) 
    8    !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions 
    9    !!            7.0  ! 1996-01  (G. Madec)  s-coordinates 
    10    !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc 
    11    !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries 
    12    !!            8.2  ! 2000-03  (G. Madec)  no slip accurate 
    13    !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90 
     6   !! History :  1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90 
    147   !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries 
    158   !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90 
     
    1710   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1811   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here 
    19    !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here 
    20    !!            3.7  ! 2015-10  (G. Madec) remove cross-land advection 
     12   !!            3.7  ! 2014-01  (G. Madec) suppression of velocity curl from in-core memory 
     13   !!             -   ! 2014-12  (G. Madec) suppression of cross land advection option 
     14   !!             -   ! 2015-10  (G. Madec) add velocity and rnf flag in argument of div_hor 
    2115   !!---------------------------------------------------------------------- 
    2216 
    2317   !!---------------------------------------------------------------------- 
    24    !!   div_cur    : Compute the horizontal divergence and relative 
    25    !!                vorticity fields 
     18   !!   div_hor    : Compute the horizontal divergence field 
    2619   !!---------------------------------------------------------------------- 
    2720   USE oce             ! ocean dynamics and tracers 
     
    2922   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean 
    3023   USE sbcrnf          ! river runoff  
    31    USE sbcisf          ! ice shelf  
     24   USE sbcisf          ! ice shelf 
     25   ! 
    3226   USE in_out_manager  ! I/O manager 
    3327   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    3933   PRIVATE 
    4034 
    41    PUBLIC   div_cur    ! routine called by step.F90 and istate.F90 
     35   PUBLIC   div_hor    ! routine called by step.F90 and istate.F90 
    4236 
    4337   !! * Substitutions 
     
    4539#  include "vectopt_loop_substitute.h90" 
    4640   !!---------------------------------------------------------------------- 
    47    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     41   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    4842   !! $Id$  
    4943   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5145CONTAINS 
    5246 
    53 #if defined key_noslip_accurate 
    54    !!---------------------------------------------------------------------- 
    55    !!   'key_noslip_accurate'   2nd order interior + 4th order at the coast 
    56    !!---------------------------------------------------------------------- 
    57  
    58    SUBROUTINE div_cur( kt ) 
     47   SUBROUTINE div_hor( kt ) 
    5948      !!---------------------------------------------------------------------- 
    60       !!                  ***  ROUTINE div_cur  *** 
     49      !!                  ***  ROUTINE div_hor  *** 
     50      !!                     
     51      !! ** Purpose :   compute the horizontal divergence at now time-step 
    6152      !! 
    62       !! ** Purpose :   compute the horizontal divergence and the relative 
    63       !!              vorticity at before and now time-step 
     53      !! ** Method  :   the now divergence is computed as : 
     54      !!         hdivn = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
     55      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)  
    6456      !! 
    65       !! ** Method  : I.  divergence : 
    66       !!         - save the divergence computed at the previous time-step 
    67       !!      (note that the Asselin filter has not been applied on hdivb) 
    68       !!         - compute the now divergence given by : 
    69       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    70       !!      correct hdiv with runoff inflow (div_rnf) and ice shelf melting (div_isf) 
    71       !!              II. vorticity : 
    72       !!         - save the curl computed at the previous time-step 
    73       !!            rotb = rotn 
    74       !!      (note that the Asselin time filter has not been applied to rotb) 
    75       !!         - compute the now curl in tensorial formalism: 
    76       !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] ) 
    77       !!         - Coastal boundary condition: 'key_noslip_accurate' defined, 
    78       !!      the no-slip boundary condition is computed using Schchepetkin 
    79       !!      and O'Brien (1996) scheme (i.e. 4th order at the coast). 
    80       !!      For example, along east coast, the one-sided finite difference 
    81       !!      approximation used for di[v] is: 
    82       !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) ) 
    83       !! 
    84       !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence 
    85       !!              - update rotb , rotn , the before & now rel. vorticity 
    86       !!---------------------------------------------------------------------- 
    87       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    88       ! 
    89       INTEGER ::   ji, jj, jk, jl           ! dummy loop indices 
    90       INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer 
    91       REAL(wp) ::  zraur, zdep              ! local scalar 
    92       REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace 
    93       REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace 
    94       !!---------------------------------------------------------------------- 
    95       ! 
    96       IF( nn_timing == 1 )  CALL timing_start('div_cur') 
    97       ! 
    98       CALL wrk_alloc( jpi  , jpj+2, zwu               ) 
    99       CALL wrk_alloc( jpi+4, jpj  , zwv, kistart = -1 ) 
    100       ! 
    101       IF( kt == nit000 ) THEN 
    102          IF(lwp) WRITE(numout,*) 
    103          IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity' 
    104          IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case' 
    105       ENDIF 
    106  
    107       !                                                ! =============== 
    108       DO jk = 1, jpkm1                                 ! Horizontal slab 
    109          !                                             ! =============== 
    110          ! 
    111          hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays 
    112          rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays 
    113          ! 
    114          !                                             ! -------- 
    115          ! Horizontal divergence                       !   div 
    116          !                                             ! -------- 
    117          DO jj = 2, jpjm1 
    118             DO ji = fs_2, fs_jpim1   ! vector opt. 
    119                hdivn(ji,jj,jk) =   & 
    120                   (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       & 
    121                    + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    & 
    122                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    123             END DO 
    124          END DO 
    125  
    126          IF( .NOT. AGRIF_Root() ) THEN 
    127             IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    128             IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west 
    129             IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north 
    130             IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south 
    131          ENDIF 
    132  
    133          !                                             ! -------- 
    134          ! relative vorticity                          !   rot  
    135          !                                             ! -------- 
    136          ! contravariant velocity (extended for lateral b.c.) 
    137          ! inside the model domain 
    138          DO jj = 1, jpj 
    139             DO ji = 1, jpi 
    140                zwu(ji,jj) = e1u(ji,jj) * un(ji,jj,jk) 
    141                zwv(ji,jj) = e2v(ji,jj) * vn(ji,jj,jk) 
    142             END DO   
    143          END DO   
    144   
    145          ! East-West boundary conditions 
    146          IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN 
    147             zwv(  0  ,:) = zwv(jpi-2,:) 
    148             zwv( -1  ,:) = zwv(jpi-3,:) 
    149             zwv(jpi+1,:) = zwv(  3  ,:) 
    150             zwv(jpi+2,:) = zwv(  4  ,:) 
    151          ELSE 
    152             zwv(  0  ,:) = 0.e0 
    153             zwv( -1  ,:) = 0.e0 
    154             zwv(jpi+1,:) = 0.e0 
    155             zwv(jpi+2,:) = 0.e0 
    156          ENDIF 
    157  
    158          ! North-South boundary conditions 
    159          IF( nperio == 3 .OR. nperio == 4 ) THEN 
    160             ! north fold ( Grid defined with a T-point pivot) ORCA 2 degre 
    161             zwu(jpi,jpj+1) = 0.e0 
    162             zwu(jpi,jpj+2) = 0.e0 
    163             DO ji = 1, jpi-1 
    164                iju = jpi - ji + 1 
    165                zwu(ji,jpj+1) = - zwu(iju,jpj-3) 
    166                zwu(ji,jpj+2) = - zwu(iju,jpj-4) 
    167             END DO 
    168          ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN 
    169             ! north fold ( Grid defined with a F-point pivot) ORCA 0.5 degre\ 
    170             zwu(jpi,jpj+1) = 0.e0 
    171             zwu(jpi,jpj+2) = 0.e0 
    172             DO ji = 1, jpi-1 
    173                iju = jpi - ji 
    174                zwu(ji,jpj  ) = - zwu(iju,jpj-1) 
    175                zwu(ji,jpj+1) = - zwu(iju,jpj-2) 
    176                zwu(ji,jpj+2) = - zwu(iju,jpj-3) 
    177             END DO 
    178             DO ji = -1, jpi+2 
    179                ijt = jpi - ji + 1 
    180                zwv(ji,jpj) = - zwv(ijt,jpj-2) 
    181             END DO 
    182             DO ji = jpi/2+1, jpi+2 
    183                ijt = jpi - ji + 1 
    184                zwv(ji,jpjm1) = - zwv(ijt,jpjm1) 
    185             END DO 
    186          ELSE 
    187             ! closed 
    188             zwu(:,jpj+1) = 0.e0 
    189             zwu(:,jpj+2) = 0.e0 
    190          ENDIF 
    191  
    192          ! relative vorticity (vertical component of the velocity curl)  
    193          DO jj = 1, jpjm1 
    194             DO ji = 1, fs_jpim1   ! vector opt. 
    195                rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      & 
    196                   &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) * r1_e1e2f(ji,jj) 
    197             END DO 
    198          END DO 
    199  
    200          ! second order accurate scheme along straight coast 
    201          DO jl = 1, npcoa(1,jk) 
    202             ii = nicoa(jl,1,jk) 
    203             ij = njcoa(jl,1,jk) 
    204             rotn(ii,ij,jk) = r1_e1e2f(ji,jj) * ( + 4. * zwv(ii+1,ij) - zwv(ii+2,ij) + 0.2 * zwv(ii+3,ij) ) 
    205          END DO 
    206          DO jl = 1, npcoa(2,jk) 
    207             ii = nicoa(jl,2,jk) 
    208             ij = njcoa(jl,2,jk) 
    209             rotn(ii,ij,jk) = r1_e1e2f(ji,jj) * (-4.*zwv(ii,ij)+zwv(ii-1,ij)-0.2*zwv(ii-2,ij)) 
    210          END DO 
    211          DO jl = 1, npcoa(3,jk) 
    212             ii = nicoa(jl,3,jk) 
    213             ij = njcoa(jl,3,jk) 
    214             rotn(ii,ij,jk) = -r1_e1e2f(ji,jj) * ( +4. * zwu(ii,ij+1) - zwu(ii,ij+2) + 0.2 * zwu(ii,ij+3) ) 
    215          END DO 
    216          DO jl = 1, npcoa(4,jk) 
    217             ii = nicoa(jl,4,jk) 
    218             ij = njcoa(jl,4,jk) 
    219             rotn(ii,ij,jk) = -r1_e1e2f(ji,jj) * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) ) 
    220          END DO 
    221          !                                             ! =============== 
    222       END DO                                           !   End of slab 
    223       !                                                ! =============== 
    224  
    225       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field) 
    226       IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
    227        
    228       ! 4. Lateral boundary conditions on hdivn and rotn 
    229       ! ---------------------------------=======---====== 
    230       CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change) 
    231       ! 
    232       CALL wrk_dealloc( jpi  , jpj+2, zwu               ) 
    233       CALL wrk_dealloc( jpi+4, jpj  , zwv, kistart = -1 ) 
    234       ! 
    235       IF( nn_timing == 1 )  CALL timing_stop('div_cur') 
    236       ! 
    237    END SUBROUTINE div_cur 
    238     
    239 #else 
    240    !!---------------------------------------------------------------------- 
    241    !!   Default option                           2nd order centered schemes 
    242    !!---------------------------------------------------------------------- 
    243  
    244    SUBROUTINE div_cur( kt ) 
    245       !!---------------------------------------------------------------------- 
    246       !!                  ***  ROUTINE div_cur  *** 
    247       !!                     
    248       !! ** Purpose :   compute the horizontal divergence and the relative 
    249       !!      vorticity at before and now time-step 
    250       !! 
    251       !! ** Method  : - Divergence: 
    252       !!      - save the divergence computed at the previous time-step 
    253       !!      (note that the Asselin filter has not been applied on hdivb) 
    254       !!      - compute the now divergence given by : 
    255       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    256       !!      correct hdiv with runoff inflow (div_rnf) 
    257       !!              - Relavtive Vorticity : 
    258       !!      - save the curl computed at the previous time-step (rotb = rotn) 
    259       !!      (note that the Asselin time filter has not been applied to rotb) 
    260       !!      - compute the now curl in tensorial formalism: 
    261       !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] ) 
    262       !!      Note: Coastal boundary condition: lateral friction set through 
    263       !!      the value of fmask along the coast (see dommsk.F90) and shlat 
    264       !!      (namelist parameter) 
    265       !! 
    266       !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence 
    267       !!              - update rotb , rotn , the before & now rel. vorticity 
     57      !! ** Action  : - update hdivn, the now horizontal divergence 
    26858      !!---------------------------------------------------------------------- 
    26959      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     
    27363      !!---------------------------------------------------------------------- 
    27464      ! 
    275       IF( nn_timing == 1 )  CALL timing_start('div_cur') 
     65      IF( nn_timing == 1 )   CALL timing_start('div_hor') 
    27666      ! 
    27767      IF( kt == nit000 ) THEN 
    27868         IF(lwp) WRITE(numout,*) 
    279          IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and' 
    280          IF(lwp) WRITE(numout,*) '~~~~~~~   relative vorticity' 
     69         IF(lwp) WRITE(numout,*) 'div_hor : horizontal velocity divergence ' 
     70         IF(lwp) WRITE(numout,*) '~~~~~~~   ' 
    28171      ENDIF 
    282  
    283       !                                                ! =============== 
    284       DO jk = 1, jpkm1                                 ! Horizontal slab 
    285          !                                             ! =============== 
    286          ! 
    287          hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays 
    288          rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays 
    289          ! 
    290          !                                             ! -------- 
    291          ! Horizontal divergence                       !   div  
    292          !                                             ! -------- 
     72      ! 
     73      DO jk = 1, jpkm1                                      !==  Horizontal divergence  ==! 
    29374         DO jj = 2, jpjm1 
    29475            DO ji = fs_2, fs_jpim1   ! vector opt. 
    295                hdivn(ji,jj,jk) =   & 
    296                   (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       & 
    297                    + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    & 
    298                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     76               hdivn(ji,jj,jk) = (  e2u(ji  ,jj) * fse3u_n(ji  ,jj,jk) * un(ji  ,jj,jk)        & 
     77                  &               - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * un(ji-1,jj,jk)        & 
     78                  &               + e1v(ji,jj  ) * fse3v_n(ji,jj  ,jk) * vn(ji,jj  ,jk)        & 
     79                  &               - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk)   )    & 
     80                  &            / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    29981            END DO   
    30082         END DO   
    301  
    30283         IF( .NOT. AGRIF_Root() ) THEN 
    303             IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    304             IF ((nbondi == -1).OR.(nbondi == 2)) hdivn(2      , :     ,jk) = 0.e0      ! west 
    305             IF ((nbondj ==  1).OR.(nbondj == 2)) hdivn(:      ,nlcj-1 ,jk) = 0.e0      ! north 
    306             IF ((nbondj == -1).OR.(nbondj == 2)) hdivn(:      ,2      ,jk) = 0.e0      ! south 
     84            IF( nbondi ==  1 .OR. nbondi == 2 )   hdivn(nlci-1,   :  ,jk) = 0._wp      ! east 
     85            IF( nbondi == -1 .OR. nbondi == 2 )   hdivn(  2   ,   :  ,jk) = 0._wp      ! west 
     86            IF( nbondj ==  1 .OR. nbondj == 2 )   hdivn(  :   ,nlcj-1,jk) = 0._wp      ! north 
     87            IF( nbondj == -1 .OR. nbondj == 2 )   hdivn(  :   ,  2   ,jk) = 0._wp      ! south 
    30788         ENDIF 
    308  
    309          !                                             ! -------- 
    310          ! relative vorticity                          !   rot  
    311          !                                             ! -------- 
    312          DO jj = 1, jpjm1 
    313             DO ji = 1, fs_jpim1   ! vector opt. 
    314                rotn(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    315                   &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
    316                   &           * fmask(ji,jj,jk) * r1_e1e2f(ji,jj) 
    317             END DO 
    318          END DO 
    319          !                                             ! =============== 
    320       END DO                                           !   End of slab 
    321       !                                                ! =============== 
    322  
    323       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field) 
    324       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
     89      END DO 
    32590      ! 
    326       CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change) 
     91      IF( ln_rnf                     )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field) 
    32792      ! 
    328       IF( nn_timing == 1 )  CALL timing_stop('div_cur') 
     93      IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    32994      ! 
    330    END SUBROUTINE div_cur 
     95      CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change) 
     96      ! 
     97      IF( nn_timing == 1 )  CALL timing_stop('div_hor') 
     98      ! 
     99   END SUBROUTINE div_hor 
    331100    
    332 #endif 
    333101   !!====================================================================== 
    334 END MODULE divcur 
     102END MODULE divhor 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r5322 r5777  
    7676      CASE ( 3 )    
    7777                      CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme 
    78       ! 
    79       CASE (-1 )                                            ! esopa: test all possibility with control print 
    80                       CALL dyn_keg     ( kt, nn_dynkeg ) 
    81                       CALL dyn_zad     ( kt ) 
    82                       CALL dyn_adv_cen2( kt ) 
    83                       CALL dyn_adv_ubs ( kt ) 
    8478      END SELECT 
    8579      ! 
     
    126120      IF( ln_dynadv_cen2 )   ioptio = ioptio + 1 
    127121      IF( ln_dynadv_ubs  )   ioptio = ioptio + 1 
    128       IF( lk_esopa       )   ioptio =          1 
    129122 
    130123      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 
     
    139132      IF( ln_dynadv_cen2 )   nadv =  2 
    140133      IF( ln_dynadv_ubs  )   nadv =  3 
    141       IF( lk_esopa       )   nadv = -1 
    142134 
    143135      IF(lwp) THEN                    ! Print the choice 
     
    151143         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used' 
    152144         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used' 
    153          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation' 
    154145      ENDIF 
    155146      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r5758 r5777  
    44   !! Ocean physics:  lateral diffusivity trends  
    55   !!===================================================================== 
    6    !! History :  9.0  !  05-11  (G. Madec)  Original code (new step architecture) 
     6   !! History :  2.0  ! 2005-11  (G. Madec)  Original code (new step architecture) 
     7   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
     8   !!                 !                                  add velocity dependent coefficient and optional read in file 
    79   !!---------------------------------------------------------------------- 
    810 
     
    1416   USE dom_oce        ! ocean space and time domain 
    1517   USE phycst         ! physical constants 
    16    USE ldfdyn_oce     ! ocean dynamics lateral physics 
    17    USE ldftra         ! ocean tracers  lateral physics 
    18    USE ldfslp         ! lateral mixing: slopes of mixing orientation 
    19    USE dynldf_bilapg  ! lateral mixing            (dyn_ldf_bilapg routine) 
    20    USE dynldf_bilap   ! lateral mixing            (dyn_ldf_bilap  routine) 
    21    USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine) 
    22    USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine) 
     18   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
     19   USE ldfslp         ! lateral diffusion: slopes of mixing orientation 
     20   USE dynldf_lap_blp ! lateral mixing   (dyn_ldf_lap & dyn_ldf_blp routines) 
     21   USE dynldf_iso     ! lateral mixing                 (dyn_ldf_iso routine ) 
    2322   USE trd_oce        ! trends: ocean variables 
    24    USE trddyn         ! trend manager: dynamics   (trd_dyn        routine) 
     23   USE trddyn         ! trend manager: dynamics   (trd_dyn      routine) 
    2524   ! 
    2625   USE prtctl         ! Print control 
     
    2827   USE lib_mpp        ! distribued memory computing library 
    2928   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    30    USE wrk_nemo        ! Memory Allocation 
    31    USE timing          ! Timing 
     29   USE wrk_nemo       ! Memory Allocation 
     30   USE timing         ! Timing 
    3231 
    3332   IMPLICIT NONE 
     
    3736   PUBLIC   dyn_ldf_init  ! called by opa  module  
    3837 
    39    INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals) 
     38   !                      ! Flag to control the type of lateral viscous operator 
     39   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in setting the operator 
     40   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral viscous trend) 
     41   !                          !!      laplacian     !    bilaplacian    ! 
     42   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator 
     43   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       ! iso-neutral or geopotential operator 
     44 
     45   INTEGER ::   nldf   ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
    4046 
    4147   !! * Substitutions 
     
    4349#  include "vectopt_loop_substitute.h90" 
    4450   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     51   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    4652   !! $Id$ 
    4753   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    6268      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf') 
    6369      ! 
    64       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
    65          CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
     70      IF( l_trddyn )   THEN                      ! temporary save of momentum trends 
     71         CALL wrk_alloc( jpi,jpj,jpk,  ztrdu, ztrdv ) 
    6672         ztrdu(:,:,:) = ua(:,:,:)  
    6773         ztrdv(:,:,:) = va(:,:,:)  
     
    7076      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
    7177      ! 
    72       CASE ( 0 )    ;   CALL dyn_ldf_lap    ( kt )      ! iso-level laplacian 
    73       CASE ( 1 )    ;   CALL dyn_ldf_iso    ( kt )      ! rotated laplacian (except dk[ dk[.] ] part) 
    74       CASE ( 2 )    ;   CALL dyn_ldf_bilap  ( kt )      ! iso-level bilaplacian 
    75 !!gm     CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian 
    76       CASE ( 4 )                                        ! iso-level laplacian + bilaplacian 
    77          CALL dyn_ldf_lap    ( kt ) 
    78          CALL dyn_ldf_bilap  ( kt ) 
    79       CASE ( 5 )                                        ! rotated laplacian + bilaplacian (s-coord) 
    80          CALL dyn_ldf_iso    ( kt ) 
    81 !!gm         CALL dyn_ldf_bilapg ( kt ) 
     78      CASE ( np_lap   )    ;   CALL dyn_ldf_lap  ( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
     79      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso  ( kt )                         ! rotated      laplacian 
     80      CASE ( np_blp   )    ;   CALL dyn_ldf_blp  ( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
    8281      ! 
    83       CASE ( -1 )                                       ! esopa: test all possibility with control print 
    84                         CALL dyn_ldf_lap    ( kt ) 
    85                         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask,   & 
    86             &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    87                         CALL dyn_ldf_iso    ( kt ) 
    88                         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask,   & 
    89             &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    90                         CALL dyn_ldf_bilap  ( kt ) 
    91                         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask,   & 
    92             &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    93 !!gm                        CALL dyn_ldf_bilapg ( kt ) 
    94                         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask,   & 
    95             &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96       ! 
    97       CASE ( -2 )                                       ! neither laplacian nor bilaplacian schemes used 
    98          IF( kt == nit000 ) THEN 
    99             IF(lwp) WRITE(numout,*) 
    100             IF(lwp) WRITE(numout,*) 'dyn_ldf : no lateral diffusion on momentum setup' 
    101             IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    102          ENDIF 
    10382      END SELECT 
    10483 
     
    10786         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    10887         CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 
    109          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
     88         CALL wrk_dealloc( jpi,jpj,jpk,  ztrdu, ztrdv ) 
    11089      ENDIF 
    11190      !                                          ! print sum trends (used for debugging) 
     
    126105      INTEGER ::   ioptio, ierr         ! temporary integers  
    127106      !!---------------------------------------------------------------------- 
    128      
     107      ! 
    129108      !                                   ! Namelist nam_dynldf: already read in ldfdyn module 
    130  
     109      ! 
    131110      IF(lwp) THEN                        ! Namelist print 
    132111         WRITE(numout,*) 
     
    134113         WRITE(numout,*) '~~~~~~~~~~~' 
    135114         WRITE(numout,*) '       Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 
    136          WRITE(numout,*) '          laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap 
    137          WRITE(numout,*) '          bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap 
    138          WRITE(numout,*) '          iso-level                   ln_dynldf_level = ', ln_dynldf_level 
    139          WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor 
    140          WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso 
     115         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap 
     116         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp 
     117         WRITE(numout,*) '          iso-level                   ln_dynldf_lev = ', ln_dynldf_lev 
     118         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor 
     119         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso 
    141120      ENDIF 
    142  
    143       !                                   ! control the consistency 
     121      !                                   ! use of lateral operator or not 
     122      nldf = np_ERROR 
    144123      ioptio = 0 
    145       IF( ln_dynldf_lap   )   ioptio = ioptio + 1 
    146       IF( ln_dynldf_bilap )   ioptio = ioptio + 1 
    147       IF( ioptio <  1 ) CALL ctl_warn( '          neither laplacian nor bilaplacian operator set for dynamics' ) 
    148       ioptio = 0 
    149       IF( ln_dynldf_level )   ioptio = ioptio + 1 
    150       IF( ln_dynldf_hor   )   ioptio = ioptio + 1 
    151       IF( ln_dynldf_iso   )   ioptio = ioptio + 1 
    152       IF( ioptio >  1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' ) 
    153  
    154       IF( ln_dynldf_iso .AND. ln_traldf_hor ) CALL ctl_stop & 
    155       &   ( 'Not sensible to use geopotential diffusion for tracers with isoneutral diffusion for dynamics' ) 
    156  
    157       !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 
    158       ierr = 0 
    159       IF ( ln_dynldf_lap ) THEN      ! laplacian operator 
    160          IF ( ln_zco ) THEN                ! z-coordinate 
    161             IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation) 
    162             IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation) 
    163             IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
     124      IF( ln_dynldf_lap )   ioptio = ioptio + 1 
     125      IF( ln_dynldf_blp )   ioptio = ioptio + 1 
     126      IF( ioptio >  1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on momentum' ) 
     127      IF( ioptio == 0   )   nldf = np_no_ldf     ! No lateral mixing operator 
     128      ! 
     129      IF( nldf /= np_no_ldf ) THEN        ! direction ==>> type of operator   
     130         ioptio = 0 
     131         IF( ln_dynldf_lev )   ioptio = ioptio + 1 
     132         IF( ln_dynldf_hor )   ioptio = ioptio + 1 
     133         IF( ln_dynldf_iso )   ioptio = ioptio + 1 
     134         IF( ioptio >  1   )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' ) 
     135         IF( ioptio == 0   )   CALL ctl_stop( '          use at least ONE direction (level/hor/iso)' ) 
     136         ! 
     137         !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 
     138         ierr = 0 
     139         IF ( ln_dynldf_lap ) THEN      ! laplacian operator 
     140            IF ( ln_zco ) THEN                ! z-coordinate 
     141               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     142               IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     143               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
     144            ENDIF 
     145            IF ( ln_zps ) THEN             ! z-coordinate with partial step 
     146               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level              (no rotation) 
     147               IF ( ln_dynldf_hor )   nldf = np_lap     ! iso-level              (no rotation) 
     148               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
     149            ENDIF 
     150            IF ( ln_sco ) THEN             ! s-coordinate 
     151               IF ( ln_dynldf_lev )   nldf = np_lap     ! iso-level = horizontal (no rotation) 
     152               IF ( ln_dynldf_hor )   nldf = np_lap_i   ! horizontal             (   rotation) 
     153               IF ( ln_dynldf_iso )   nldf = np_lap_i   ! iso-neutral            (   rotation) 
     154            ENDIF 
    164155         ENDIF 
    165          IF ( ln_zps ) THEN             ! z-coordinate 
    166             IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed 
    167             IF ( ln_dynldf_hor   )   nldf = 0      ! horizontal (no rotation) 
    168             IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
     156         ! 
     157         IF( ln_dynldf_blp ) THEN          ! bilaplacian operator 
     158            IF ( ln_zco ) THEN                ! z-coordinate 
     159               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     160               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level = horizontal (no rotation) 
     161               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
     162            ENDIF 
     163            IF ( ln_zps ) THEN             ! z-coordinate with partial step 
     164               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation) 
     165               IF ( ln_dynldf_hor )   nldf = np_blp     ! iso-level              (no rotation) 
     166               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
     167            ENDIF 
     168            IF ( ln_sco ) THEN             ! s-coordinate 
     169               IF ( ln_dynldf_lev )   nldf = np_blp     ! iso-level              (no rotation) 
     170               IF ( ln_dynldf_hor )   ierr = 2          ! horizontal             (   rotation) 
     171               IF ( ln_dynldf_iso )   ierr = 2          ! iso-neutral            (   rotation) 
     172            ENDIF 
    169173         ENDIF 
    170          IF ( ln_sco ) THEN             ! s-coordinate 
    171             IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation) 
    172             IF ( ln_dynldf_hor   )   nldf = 1      ! horizontal (   rotation) 
    173             IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation) 
    174          ENDIF 
     174         ! 
     175         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 
     176         ! 
     177         IF( nldf == np_lap_i )   l_ldfslp = .TRUE.      ! rotation require the computation of the slopes 
     178         ! 
    175179      ENDIF 
    176  
    177       IF( ln_dynldf_bilap ) THEN      ! bilaplacian operator 
    178          IF ( ln_zco ) THEN                ! z-coordinate 
    179             IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation) 
    180             IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation) 
    181             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    182          ENDIF 
    183          IF ( ln_zps ) THEN             ! z-coordinate 
    184             IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed  
    185             IF ( ln_dynldf_hor   )   nldf = 2      ! horizontal (no rotation) 
    186             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    187          ENDIF 
    188          IF ( ln_sco ) THEN             ! s-coordinate 
    189             IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation) 
    190             IF ( ln_dynldf_hor   )   nldf = 3      ! horizontal (   rotation) 
    191             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    192          ENDIF 
    193       ENDIF 
    194        
    195       IF( ln_dynldf_lap .AND. ln_dynldf_bilap ) THEN  ! mixed laplacian and bilaplacian operators 
    196          IF ( ln_zco ) THEN                ! z-coordinate 
    197             IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation) 
    198             IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation) 
    199             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    200          ENDIF 
    201          IF ( ln_zps ) THEN             ! z-coordinate 
    202             IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed  
    203             IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation) 
    204             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    205          ENDIF 
    206          IF ( ln_sco ) THEN             ! s-coordinate 
    207             IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation) 
    208             IF ( ln_dynldf_hor   )   nldf = 5      ! horizontal (   rotation) 
    209             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation) 
    210          ENDIF 
    211       ENDIF 
    212  
    213       IF( lk_esopa )                 nldf = -1     ! esopa test 
    214  
    215       IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' ) 
    216       IF( ierr == 2 )   CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' ) 
    217       IF( nldf == 1 .OR. nldf == 3 )   l_ldfslp = .TRUE.    ! the rotation needs slope computation 
    218180 
    219181      IF(lwp) THEN 
    220182         WRITE(numout,*) 
    221          IF( nldf == -2 )   WRITE(numout,*) '              neither laplacian nor bilaplacian schemes used' 
    222          IF( nldf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used' 
    223          IF( nldf ==  0 )   WRITE(numout,*) '              laplacian operator' 
    224          IF( nldf ==  1 )   WRITE(numout,*) '              rotated laplacian operator' 
    225          IF( nldf ==  2 )   WRITE(numout,*) '              bilaplacian operator' 
    226          IF( nldf ==  3 )   WRITE(numout,*) '              rotated bilaplacian' 
    227          IF( nldf ==  4 )   WRITE(numout,*) '              laplacian and bilaplacian operators' 
    228          IF( nldf ==  5 )   WRITE(numout,*) '              rotated laplacian and bilaplacian operators' 
     183         IF( nldf == np_no_ldf )   WRITE(numout,*) '              NO lateral viscosity' 
     184         IF( nldf == np_lap    )   WRITE(numout,*) '              iso-level laplacian operator' 
     185         IF( nldf == np_lap_i  )   WRITE(numout,*) '              rotated laplacian operator with iso-level background' 
     186         IF( nldf == np_blp    )   WRITE(numout,*) '              iso-level bi-laplacian operator' 
    229187      ENDIF 
    230188      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r5758 r5777  
    22   !!====================================================================== 
    33   !!                     ***  MODULE  dynldf_iso  *** 
    4    !! Ocean dynamics:  lateral viscosity trend 
     4   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator) 
    55   !!====================================================================== 
    66   !! History :  OPA  !  97-07  (G. Madec)  Original code 
     
    88   !!             -   !  2004-08  (C. Talandier) New trends organization 
    99   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion 
     10   !!            3.7  !  2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
     11   !!                 !                                   add velocity dependent coefficient and optional read in file 
    1012   !!---------------------------------------------------------------------- 
    1113 
     
    1719   USE oce             ! ocean dynamics and tracers 
    1820   USE dom_oce         ! ocean space and time domain 
    19    USE ldfdyn_oce      ! ocean dynamics lateral physics 
    20    USE ldftra          ! lateral physics: eddy diffusivity & EIV coefficients 
     21   USE ldfdyn          ! lateral diffusion: eddy viscosity coef. 
     22   USE ldftra          ! lateral physics: eddy diffusivity 
    2123   USE zdf_oce         ! ocean vertical physics 
    2224   USE ldfslp          ! iso-neutral slopes  
    2325   ! 
    24    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2526   USE in_out_manager  ! I/O manager 
    2627   USE lib_mpp         ! MPP library 
     28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2729   USE prtctl          ! Print control 
    2830   USE wrk_nemo        ! Memory Allocation 
     
    4042   !! * Substitutions 
    4143#  include "domzgr_substitute.h90" 
    42 #  include "ldfdyn_substitute.h90" 
    4344#  include "vectopt_loop_substitute.h90" 
    4445   !!---------------------------------------------------------------------- 
     
    8182      !!      horizontal fluxes associated with the rotated lateral mixing: 
    8283      !!      u-component: 
    83       !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ] 
    84       !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 
    85       !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ] 
    86       !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ] 
     84      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ] 
     85      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 
     86      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
     87      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ] 
    8788      !!      v-component: 
    88       !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ] 
    89       !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ] 
    90       !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ] 
    91       !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 
     89      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ] 
     90      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ] 
     91      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
     92      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 
    9293      !!      take the horizontal divergence of the fluxes: 
    9394      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  } 
     
    108109      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    109110      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars 
    110       REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      - 
     111      REAL(wp) ::   zmskt, zmskf                                     !   -      - 
    111112      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      - 
    112113      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
     
    127128      ENDIF 
    128129 
    129       ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 
     130!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED 
     131      ! s-coordinate: Iso-level diffusion on momentum but not on tracer 
    130132      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 
    131133         ! 
     
    133135            DO jj = 2, jpjm1 
    134136               DO ji = 2, jpim1 
    135                   uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
    136                   vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
    137                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
    138                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
     137                  uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
     138                  vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
     139                  wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 
     140                  wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 
    139141               END DO 
    140142            END DO 
     
    181183            DO jj = 2, jpjm1 
    182184               DO ji = fs_2, jpi   ! vector opt. 
    183                   zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 
    184  
    185                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    186                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
     185                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj) 
     186 
     187                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     & 
     188                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
    187189 
    188190                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    189191    
     192                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    & 
     193                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      & 
     194                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     195               END DO 
     196            END DO 
     197         ELSE                   ! other coordinate system (zco or sco) : e3t 
     198            DO jj = 2, jpjm1 
     199               DO ji = fs_2, jpi   ! vector opt. 
     200                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj) 
     201 
     202                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     & 
     203                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp ) 
     204 
     205                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     206 
    190207                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    191208                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
     
    193210               END DO 
    194211            END DO 
    195          ELSE                   ! other coordinate system (zco or sco) : e3t 
    196             DO jj = 2, jpjm1 
    197                DO ji = fs_2, jpi   ! vector opt. 
    198                   zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
    199  
    200                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    201                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    202  
    203                   zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    204  
    205                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    206                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
    207                      &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
    208                END DO 
    209             END DO 
    210212         ENDIF 
    211213 
     
    213215         DO jj = 1, jpjm1 
    214216            DO ji = 1, fs_jpim1   ! vector opt. 
    215                zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
    216  
    217                zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
    218                   &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
     217               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj) 
     218 
     219               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     & 
     220                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
    219221 
    220222               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     
    234236         DO jj = 2, jpjm1 
    235237            DO ji = 1, fs_jpim1   ! vector opt. 
    236                zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
    237  
    238                zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
    239                   &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
     238               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj) 
     239 
     240               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     & 
     241                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
    240242 
    241243               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    242244 
    243                zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
    244                   &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     & 
    245                   &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
     245               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    & 
     246                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      & 
     247                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
    246248            END DO 
    247249         END DO 
     
    251253            DO jj = 2, jpj 
    252254               DO ji = 1, fs_jpim1   ! vector opt. 
    253                   zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 
     255                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj) 
     256 
     257                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     & 
     258                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
     259 
     260                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     261 
     262                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    & 
     263                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      & 
     264                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     265               END DO 
     266            END DO 
     267         ELSE                   ! other coordinate system (zco or sco) : e3t 
     268            DO jj = 2, jpj 
     269               DO ji = 1, fs_jpim1   ! vector opt. 
     270                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj) 
    254271 
    255272                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     
    263280               END DO 
    264281            END DO 
    265          ELSE                   ! other coordinate system (zco or sco) : e3t 
    266             DO jj = 2, jpj 
    267                DO ji = 1, fs_jpim1   ! vector opt. 
    268                   zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
    269  
    270                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    271                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    272  
    273                   zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    274  
    275                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
    276                      &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
    277                      &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
    278                END DO 
    279             END DO 
    280282         ENDIF 
    281283 
     
    283285         ! Second derivative (divergence) and add to the general trend 
    284286         ! ----------------------------------------------------------- 
    285  
    286287         DO jj = 2, jpjm1 
    287             DO ji = 2, jpim1          !! Question vectop possible??? !!bug 
    288                ! volume elements 
    289                zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    290                zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    291                ! horizontal component of isopycnal momentum diffusive trends 
    292                zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   & 
    293                   &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu 
    294                zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   & 
    295                   &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv 
    296                ! add the trends to the general trends 
    297                ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 
    298                va (ji,jj,jk) = va (ji,jj,jk) + zvah 
     288            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug 
     289               ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj  )    & 
     290                  &                          + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     291               va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj  ) - zivf(ji-1,jj)    & 
     292                  &                          + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    299293            END DO 
    300294         END DO 
     
    362356               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
    363357                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
    364                zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   & 
    365                              + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. ) 
     358               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   & 
     359                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. ) 
    366360 
    367361               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 
     
    409403         DO jk = 1, jpkm1 
    410404            DO ji = 2, jpim1 
    411                ! volume elements 
    412                zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 
    413                zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 
    414                ! part of the k-component of isopycnal momentum diffusive trends 
    415                zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 
    416                zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 
    417                ! add the trends to the general trends 
    418                ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 
    419                va(ji,jj,jk) = va(ji,jj,jk) + zvav 
     405               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     406               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    420407            END DO 
    421408         END DO 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90

    r5775 r5777  
    1 MODULE dynldf_lap 
     1MODULE dynldf_lap_blp 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  dynldf_lap  *** 
    4    !! Ocean dynamics:  lateral viscosity trend 
     3   !!                   ***  MODULE  dynldf_lap_blp  *** 
     4   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian) 
    55   !!====================================================================== 
    66   !! History :  OPA  ! 1990-09 (G. Madec) Original code 
     
    99   !!   NEMO     1.0  ! 2002-06 (G. Madec)  F90: Free form and module 
    1010   !!             -   ! 2004-08 (C. Talandier) New trends organization 
     11   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
     12   !!                 !                                  add velocity dependent coefficient and optional read in file 
    1113   !!---------------------------------------------------------------------- 
    1214 
    1315   !!---------------------------------------------------------------------- 
    14    !!   dyn_ldf_lap  : update the momentum trend with the lateral diffusion 
    15    !!                  using an iso-level harmonic operator 
     16   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator 
     17   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator 
    1618   !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and tracers 
    18    USE dom_oce         ! ocean space and time domain 
    19    USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    20    USE zdf_oce         ! ocean vertical physics 
     19   USE oce            ! ocean dynamics and tracers 
     20   USE dom_oce        ! ocean space and time domain 
     21   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
     22   USE ldfslp         ! iso-neutral slopes  
     23   USE zdf_oce        ! ocean vertical physics 
    2124   ! 
    22    USE in_out_manager  ! I/O manager 
    23    USE timing          ! Timing 
     25   USE in_out_manager ! I/O manager 
     26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     27   USE wrk_nemo       ! Memory Allocation 
     28   USE timing         ! Timing 
    2429 
    2530   IMPLICIT NONE 
    2631   PRIVATE 
    2732 
    28    PUBLIC dyn_ldf_lap  ! called by step.F90 
     33   PUBLIC dyn_ldf_lap  ! called by dynldf.F90 
     34   PUBLIC dyn_ldf_blp  ! called by dynldf.F90 
    2935 
    3036   !! * Substitutions 
    3137#  include "domzgr_substitute.h90" 
    32 #  include "ldfdyn_substitute.h90" 
    3338#  include "vectopt_loop_substitute.h90" 
    3439   !!---------------------------------------------------------------------- 
    35    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     40   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    3641   !! $Id$  
    3742   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    3944CONTAINS 
    4045 
    41    SUBROUTINE dyn_ldf_lap( kt ) 
     46   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) 
    4247      !!---------------------------------------------------------------------- 
    4348      !!                     ***  ROUTINE dyn_ldf_lap  *** 
    4449      !!                        
    45       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    46       !!      trend and add it to the general trend of tracer equation. 
     50      !! ** Purpose :   Compute the before horizontal momentum diffusive  
     51      !!      trend and add it to the general trend of momentum equation. 
    4752      !! 
    48       !! ** Method  :   The before horizontal momentum diffusion trend is an 
    49       !!      harmonic operator (laplacian type) which separates the divergent 
    50       !!      and rotational parts of the flow. 
    51       !!      Its horizontal components are computed as follow: 
    52       !!         difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb] 
    53       !!         difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb] 
    54       !!      in the rotational part of the diffusion. 
    55       !!      Add this before trend to the general trend (ua,va): 
    56       !!            (ua,va) = (ua,va) + (diffu,diffv) 
     53      !! ** Method  :   The Laplacian operator apply on horizontal velocity is  
     54      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
    5755      !! 
    58       !! ** Action : - Update (ua,va) with the iso-level harmonic mixing trend 
     56      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 
    5957      !!---------------------------------------------------------------------- 
    60       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     58      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     59      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     60      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
     61      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2] 
    6162      ! 
    62       INTEGER  ::   ji, jj, jk             ! dummy loop indices 
    63       REAL(wp) ::   zua, zva, ze2u, ze1v   ! local scalars 
     63      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     64      REAL(wp) ::   zsign        ! local scalars 
     65      REAL(wp) ::   zua, zva     ! local scalars 
     66      REAL(wp), POINTER, DIMENSION(:,:) ::  zcur, zdiv 
    6467      !!---------------------------------------------------------------------- 
    6568      ! 
    66       IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_lap') 
     69      IF( kt == nit000 .AND. lwp ) THEN 
     70         WRITE(numout,*) 
     71         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 
     72         WRITE(numout,*) '~~~~~~~ ' 
     73      ENDIF 
    6774      ! 
    68       IF( kt == nit000 ) THEN 
    69          IF(lwp) WRITE(numout,*) 
    70          IF(lwp) WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator' 
    71          IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
     75      IF( nn_timing == 1 )   CALL timing_start('dyn_ldf_lap') 
     76      ! 
     77      CALL wrk_alloc( jpi, jpj, zcur, zdiv )  
     78      ! 
     79      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign 
     80      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0) 
    7281      ENDIF 
     82      ! 
    7383      !                                                ! =============== 
    7484      DO jk = 1, jpkm1                                 ! Horizontal slab 
    7585         !                                             ! =============== 
    76          DO jj = 2, jpjm1 
     86         DO jj = 2, jpj 
     87            DO ji = fs_2, jpi   ! vector opt. 
     88               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
     89               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * fse3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      & 
     90                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)                & 
     91                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) * fmask(ji-1,jj-1,jk) 
     92               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
     93               zdiv(ji,jj)     = ahmt(ji,jj,jk) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )                     * tmask(ji,jj,jk)    & 
     94                  &     * (  e2u(ji,jj)*fse3u(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * pub(ji-1,jj,jk)    & 
     95                  &        + e1v(ji,jj)*fse3v(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
     96            END DO   
     97         END DO   
     98         ! 
     99         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div ) 
    77100            DO ji = fs_2, fs_jpim1   ! vector opt. 
    78                ze2u = rotb (ji,jj,jk) * fsahmf(ji,jj,jk) * fse3f(ji,jj,jk) 
    79                ze1v = hdivb(ji,jj,jk) * fsahmt(ji,jj,jk) 
    80                ! horizontal diffusive trends 
    81                zua = - ( ze2u - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk)*fse3f(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
    82                      + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) - ze1v                   ) * r1_e1u(ji,jj) 
    83  
    84                zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk)*fse3f(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
    85                      + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v                   ) * r1_e2v(ji,jj) 
    86  
    87                ! add it to the general momentum trends 
    88                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    89                va(ji,jj,jk) = va(ji,jj,jk) + zva 
     101               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                   & 
     102                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) /  ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
     103                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
     104                  ! 
     105               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                   & 
     106                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) /  ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
     107                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
    90108            END DO 
    91109         END DO 
     
    93111      END DO                                           !   End of slab 
    94112      !                                                ! =============== 
     113      CALL wrk_dealloc( jpi, jpj, zcur, zdiv )  
     114      ! 
    95115      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap') 
    96116      ! 
    97117   END SUBROUTINE dyn_ldf_lap 
    98118 
     119 
     120   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) 
     121      !!---------------------------------------------------------------------- 
     122      !!                 ***  ROUTINE dyn_ldf_blp  *** 
     123      !!                     
     124      !! ** Purpose :   Compute the before lateral momentum viscous trend  
     125      !!              and add it to the general trend of momentum equation. 
     126      !! 
     127      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian 
     128      !!      operator applied to before field (forward in time). 
     129      !!      It is computed by two successive calls to dyn_ldf_lap routine 
     130      !! 
     131      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion 
     132      !!---------------------------------------------------------------------- 
     133      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields 
     135      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend 
     136      ! 
     137      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zulap, zvlap   ! laplacian at u- and v-point 
     138      !!---------------------------------------------------------------------- 
     139      ! 
     140      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_blp') 
     141      ! 
     142      CALL wrk_alloc( jpi, jpj, jpk, zulap, zvlap )  
     143      ! 
     144      IF( kt == nit000 )  THEN 
     145         IF(lwp) WRITE(numout,*) 
     146         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' 
     147         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     148      ENDIF 
     149      ! 
     150      zulap(:,:,:) = 0._wp 
     151      zvlap(:,:,:) = 0._wp 
     152      ! 
     153      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap) 
     154      ! 
     155      CALL lbc_lnk( zulap(:,:,:) , 'U', -1. )             ! Lateral boundary conditions 
     156      CALL lbc_lnk( zvlap(:,:,:) , 'V', -1. ) 
     157      ! 
     158      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta) 
     159      ! 
     160      CALL wrk_dealloc( jpi, jpj, jpk, zulap, zvlap )  
     161      ! 
     162      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_blp') 
     163      ! 
     164   END SUBROUTINE dyn_ldf_blp 
     165 
    99166   !!====================================================================== 
    100 END MODULE dynldf_lap 
     167END MODULE dynldf_lap_blp 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5770 r5777  
    7777      !!             Note that as all external forcing a time averaging over a two rdt 
    7878      !!             period is used to prevent the divergence of odd and even time step. 
    79       !! 
    80       !! N.B. : When key_esopa is used all the scheme are tested, regardless  
    81       !!        of the physical meaning of the results.  
    8279      !!---------------------------------------------------------------------- 
    8380      ! 
     
    121118               DO ji = fs_2, fs_jpim1   ! vector opt. 
    122119                  spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    & 
    123                      &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) /e1u(ji,jj) 
     120                     &                      + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj) 
    124121                  spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    & 
    125                      &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) /e2v(ji,jj) 
     122                     &                      + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj) 
    126123               END DO 
    127124            END DO 
     
    135132            DO jj = 2, jpjm1                         ! add tide potential forcing 
    136133               DO ji = fs_2, fs_jpim1   ! vector opt. 
    137                   spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
    138                   spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     134                  spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 
     135                  spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 
    139136               END DO  
    140137            END DO 
     
    149146            DO jj = 2, jpjm1 
    150147               DO ji = fs_2, fs_jpim1   ! vector opt. 
    151                   spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    152                   spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     148                  spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     149                  spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    153150               END DO 
    154151            END DO 
     
    176173      CASE (  2 )   ;   CALL dyn_spg_flt( kt, kindic )      ! filtered 
    177174      !                                                     
    178       CASE ( -1 )                                ! esopa: test all possibility with control print 
    179                         CALL dyn_spg_exp( kt ) 
    180                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, & 
    181          &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    182                         CALL dyn_spg_ts ( kt ) 
    183                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, & 
    184          &                           tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    185                         CALL dyn_spg_flt( kt, kindic ) 
    186                         CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, & 
    187          &                            tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    188175      END SELECT 
    189176      !                     
     
    248235      IF(lk_dynspg_flt)   ioptio = ioptio + 1 
    249236      ! 
    250       IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
     237      IF(  ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    251238           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 
    252239      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   & 
    253240           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 
    254241      ! 
    255       IF( lk_esopa     )   nspg = -1 
    256242      IF( lk_dynspg_exp)   nspg =  0 
    257243      IF( lk_dynspg_ts )   nspg =  1 
    258244      IF( lk_dynspg_flt)   nspg =  2 
    259245      ! 
    260       IF( lk_esopa     )   nspg = -1 
    261       ! 
    262246      IF(lwp) THEN 
    263247         WRITE(numout,*) 
    264          IF( nspg == -1 )   WRITE(numout,*) '     ESOPA test All scheme used' 
    265248         IF( nspg ==  0 )   WRITE(numout,*) '     explicit free surface' 
    266249         IF( nspg ==  1 )   WRITE(numout,*) '     free surface with time splitting scheme' 
     
    268251      ENDIF 
    269252 
    270 #if defined key_dynspg_flt || defined key_esopa 
     253#if defined key_dynspg_flt 
    271254      CALL solver_init( nit000 )   ! Elliptic solver initialisation 
    272255#endif 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r4990 r5777  
    77   !!            3.2  !  2009-06  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 
    88   !!---------------------------------------------------------------------- 
    9 #if defined key_dynspg_exp   ||   defined key_esopa 
     9#if defined key_dynspg_exp 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_dynspg_exp'                              explicit free surface 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r5770 r5777  
    1616   !!             -   !  2014-12  (G. Madec) remove cross-land advection (cla) 
    1717   !!---------------------------------------------------------------------- 
    18 #if defined key_dynspg_flt   ||   defined key_esopa   
     18#if defined key_dynspg_flt   
    1919   !!---------------------------------------------------------------------- 
    2020   !!   'key_dynspg_flt'                              filtered free surface 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r4486 r5777  
    1717    
    1818   !                                                       !!! Surface pressure gradient logicals 
    19 #if   defined key_dynspg_exp  ||  defined key_esopa 
     19#if   defined key_dynspg_exp 
    2020   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .TRUE.   !: Explicit free surface flag 
    2121#else 
    2222   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .FALSE.  !: Explicit free surface flag 
    2323#endif 
    24 #if   defined key_dynspg_ts   ||  defined key_esopa 
     24#if   defined key_dynspg_ts 
    2525   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .TRUE.   !: Free surface with time splitting flag 
    2626#else 
    2727   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .FALSE.  !: Free surface with time splitting flag 
    2828#endif 
    29 #if   defined key_dynspg_flt  ||  defined key_esopa 
     29#if   defined key_dynspg_flt 
    3030   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .TRUE.   !: Filtered free surface cst volume flag 
    3131#else 
    3232   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .FALSE.  !: Filtered free surface cst volume flag 
    3333#endif 
    34  
    3534  !                                                                         !!! Time splitting scheme (key_dynspg_ts)  
    3635   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_e, ssha_e   ! sea surface heigth (now, after) 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5737 r5777  
    1212   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1313   !!--------------------------------------------------------------------- 
    14 #if defined key_dynspg_ts   ||   defined key_esopa 
     14#if defined key_dynspg_ts 
    1515   !!---------------------------------------------------------------------- 
    1616   !!   'key_dynspg_ts'         split explicit free surface 
     
    9898      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
    9999 
    100       IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    101                                                     &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
     100      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     101         &                          ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
    102102 
    103103      dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
     
    220220      ! 
    221221      IF ( kt == nit000 .OR. lk_vvl ) THEN 
    222          IF ( ln_dynvor_een_old ) THEN 
    223             DO jj = 1, jpjm1 
    224                DO ji = 1, jpim1 
    225                   zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
    226                         &          ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp   
    227                   IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
    228                END DO 
    229             END DO 
     222         IF ( ln_dynvor_een ) THEN              !==  EEN scheme  ==! 
     223            SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point 
     224            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     225               DO jj = 1, jpjm1 
     226                  DO ji = 1, jpim1 
     227                     zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
     228                        &             ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp   
     229                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 
     230                  END DO 
     231               END DO 
     232            CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     233               DO jj = 1, jpjm1 
     234                  DO ji = 1, jpim1 
     235                     zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     & 
     236                        &             ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   & 
     237                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
     238                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
     239                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 
     240                  END DO 
     241               END DO 
     242            END SELECT 
    230243            CALL lbc_lnk( zwz, 'F', 1._wp ) 
    231             zwz(:,:) = ff(:,:) * zwz(:,:) 
    232  
     244            ! 
    233245            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    234246            DO jj = 2, jpj 
    235                DO ji = fs_2, jpi   ! vector opt. 
     247               DO ji = 2, jpi 
    236248                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    237249                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     
    240252               END DO 
    241253            END DO 
    242          ELSE IF ( ln_dynvor_een ) THEN 
    243             DO jj = 1, jpjm1 
    244                DO ji = 1, jpim1 
    245                   zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                     & 
    246                         &          ht(ji  ,jj  ) + ht(ji+1,jj  )   )                   & 
    247                         &      / ( MAX( 1.0_wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
    248                         &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
    249                   IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
    250                END DO 
    251             END DO 
    252             CALL lbc_lnk( zwz, 'F', 1._wp ) 
    253             zwz(:,:) = ff(:,:) * zwz(:,:) 
    254  
    255             ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
    256             DO jj = 2, jpj 
    257                DO ji = fs_2, jpi   ! vector opt. 
    258                   ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
    259                   ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
    260                   ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
    261                   ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
    262                END DO 
    263             END DO 
    264          ELSE 
     254            ! 
     255         ELSE                                !== all other schemes (ENE, ENS, MIX) 
    265256            zwz(:,:) = 0._wp 
    266             zhf(:,:) = 0. 
     257            zhf(:,:) = 0._wp 
    267258            IF ( .not. ln_sco ) THEN 
     259 
     260!!gm  agree the JC comment  : this should be done in a much clear way 
     261 
    268262! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    269263!     Set it to zero for the time being  
     
    277271 
    278272            DO jj = 1, jpjm1 
    279                zhf(:,jj) = zhf(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     273               zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    280274            END DO 
    281275 
     
    298292      ! If forward start at previous time step, and centered integration,  
    299293      ! then update averaging weights: 
    300       IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 
     294      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 
    301295         ll_fw_start=.FALSE. 
    302296         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
     
    361355         END DO 
    362356         ! 
    363       ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN  ! enstrophy and energy conserving scheme 
     357      ELSEIF ( ln_dynvor_een ) THEN  ! enstrophy and energy conserving scheme 
    364358         DO jj = 2, jpjm1 
    365359            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    710704            END DO 
    711705            ! 
    712          ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==! 
     706         ELSEIF ( ln_dynvor_een ) THEN !==  energy and enstrophy conserving scheme  ==! 
    713707            DO jj = 2, jpjm1 
    714708               DO ji = fs_2, fs_jpim1   ! vector opt. 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r5737 r5777  
    1616   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    1717   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity  
     18   !!             -   ! 2014-06  (G. Madec) suppression of velocity curl from in-core memory 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    2223   !!       vor_ens  : enstrophy conserving scheme       (ln_dynvor_ens=T) 
    2324   !!       vor_ene  : energy conserving scheme          (ln_dynvor_ene=T) 
    24    !!       vor_mix  : mixed enstrophy/energy conserving (ln_dynvor_mix=T) 
    2525   !!       vor_een  : energy and enstrophy conserving   (ln_dynvor_een=T) 
    2626   !!   dyn_vor_init : set and control of the different vorticity option 
     
    3232   USE trd_oce        ! trends: ocean variables 
    3333   USE trddyn         ! trend manager: dynamics 
     34   ! 
    3435   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3536   USE prtctl         ! Print control 
     
    4445 
    4546   PUBLIC   dyn_vor        ! routine called by step.F90 
    46    PUBLIC   dyn_vor_init   ! routine called by opa.F90 
     47   PUBLIC   dyn_vor_init   ! routine called by nemogcm.F90 
    4748 
    4849   !                                   !!* Namelist namdyn_vor: vorticity term 
    49    LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme 
    50    LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme 
    51    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme 
    52    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme 
    53    LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 
    54  
    55    INTEGER ::   nvor = 0   ! type of vorticity trend used 
    56    INTEGER ::   ncor = 1   ! coriolis 
    57    INTEGER ::   nrvm = 2   ! =2 relative vorticity ; =3 metric term 
    58    INTEGER ::   ntot = 4   ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term 
    59  
     50   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: energy conserving scheme    (ENE) 
     51   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme (ENS) 
     52   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                (MIX) 
     53   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme (EEN) 
     54   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     55   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     56 
     57   INTEGER ::   nvor_scheme        ! choice of the type of advection scheme 
     58   !                               ! associated indices: 
     59   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme 
     60   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 2   ! ENS scheme 
     61   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 3   ! MIX scheme 
     62   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
     63 
     64   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     65   !                               ! associated indices: 
     66   INTEGER, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
     67   INTEGER, PARAMETER ::   np_RVO = 2         ! relative vorticity 
     68   INTEGER, PARAMETER ::   np_MET = 3         ! metric term 
     69   INTEGER, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity) 
     70   INTEGER, PARAMETER ::   np_CME = 5         ! Coriolis + metric term 
     71    
     72   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     73   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 
     74   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 
     75    
    6076   !! * Substitutions 
    6177#  include "domzgr_substitute.h90" 
    6278#  include "vectopt_loop_substitute.h90" 
    6379   !!---------------------------------------------------------------------- 
    64    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     80   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    6581   !! $Id$ 
    6682   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    87103      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    88104      ! 
    89       !                                          ! vorticity term  
    90       SELECT CASE ( nvor )                       ! compute the vorticity trend and add it to the general trend 
    91       ! 
    92       CASE ( -1 )                                      ! esopa: test all possibility with control print 
    93          CALL vor_ene( kt, ntot, ua, va ) 
    94          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 
    95             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96          CALL vor_ens( kt, ntot, ua, va ) 
    97          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 
    98             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    99          CALL vor_mix( kt ) 
    100          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 
    101             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    102          CALL vor_een( kt, ntot, ua, va ) 
    103          CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 
    104             &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    105          ! 
    106       CASE ( 0 )                                       ! energy conserving scheme 
    107          IF( l_trddyn )   THEN 
    108             ztrdu(:,:,:) = ua(:,:,:) 
    109             ztrdv(:,:,:) = va(:,:,:) 
    110             CALL vor_ene( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     105      SELECT CASE ( nvor_scheme )               !==  vorticity trend added to the general trend  ==! 
     106      ! 
     107      CASE ( np_ENE )                                 !* energy conserving scheme 
     108         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     109            ztrdu(:,:,:) = ua(:,:,:) 
     110            ztrdv(:,:,:) = va(:,:,:) 
     111            CALL vor_ene( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
    111112            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    112113            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    114115            ztrdu(:,:,:) = ua(:,:,:) 
    115116            ztrdv(:,:,:) = va(:,:,:) 
    116             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend 
     117            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend 
    117118            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    118119            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    119120            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    120121         ELSE 
    121             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity 
    122          ENDIF 
    123          ! 
    124       CASE ( 1 )                                       ! enstrophy conserving scheme 
    125          IF( l_trddyn )   THEN     
    126             ztrdu(:,:,:) = ua(:,:,:) 
    127             ztrdv(:,:,:) = va(:,:,:) 
    128             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     122            CALL vor_ene( kt, ntot, ua, va )                ! total vorticity trend 
     123         ENDIF 
     124         ! 
     125      CASE ( np_ENS )                                 !* enstrophy conserving scheme 
     126         IF( l_trddyn ) THEN                                ! trend diagnostics: splitthe trend in two     
     127            ztrdu(:,:,:) = ua(:,:,:) 
     128            ztrdv(:,:,:) = va(:,:,:) 
     129            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
    129130            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    130131            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    132133            ztrdu(:,:,:) = ua(:,:,:) 
    133134            ztrdv(:,:,:) = va(:,:,:) 
    134             CALL vor_ens( kt, ncor, ua, va )                ! planetary vorticity trend 
     135            CALL vor_ens( kt, ncor, ua, va )                      ! planetary vorticity trend 
    135136            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    136137            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    137138            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    138139         ELSE 
    139             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity 
    140          ENDIF 
    141          ! 
    142       CASE ( 2 )                                       ! mixed ene-ens scheme 
    143          IF( l_trddyn )   THEN 
    144             ztrdu(:,:,:) = ua(:,:,:) 
    145             ztrdv(:,:,:) = va(:,:,:) 
    146             CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
     140            CALL vor_ens( kt, ntot, ua, va )                ! total vorticity trend 
     141         ENDIF 
     142         ! 
     143      CASE ( np_MIX )                                 !* mixed ene-ens scheme 
     144         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     145            ztrdu(:,:,:) = ua(:,:,:) 
     146            ztrdv(:,:,:) = va(:,:,:) 
     147            CALL vor_ens( kt, nrvm, ua, va )                      ! relative vorticity or metric trend (ens) 
    147148            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    148149            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    150151            ztrdu(:,:,:) = ua(:,:,:) 
    151152            ztrdv(:,:,:) = va(:,:,:) 
    152             CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
     153            CALL vor_ene( kt, ncor, ua, va )                      ! planetary vorticity trend (ene) 
    153154            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    154155            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    155156            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    156157         ELSE 
    157             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene) 
    158          ENDIF 
    159          ! 
    160       CASE ( 3 )                                       ! energy and enstrophy conserving scheme 
    161          IF( l_trddyn )   THEN 
    162             ztrdu(:,:,:) = ua(:,:,:) 
    163             ztrdv(:,:,:) = va(:,:,:) 
    164             CALL vor_een( kt, nrvm, ua, va )                ! relative vorticity or metric trend 
     158            CALL vor_ens( kt, nrvm, ua, va )                ! relative vorticity or metric trend (ens) 
     159            CALL vor_ene( kt, ncor, ua, va )                ! planetary vorticity trend (ene) 
     160        ENDIF 
     161         ! 
     162      CASE ( np_EEN )                                 !* energy and enstrophy conserving scheme 
     163         IF( l_trddyn ) THEN                                ! trend diagnostics: split the trend in two 
     164            ztrdu(:,:,:) = ua(:,:,:) 
     165            ztrdv(:,:,:) = va(:,:,:) 
     166            CALL vor_een( kt, nrvm, ua, va )                      ! relative vorticity or metric trend 
    165167            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    166168            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     
    168170            ztrdu(:,:,:) = ua(:,:,:) 
    169171            ztrdv(:,:,:) = va(:,:,:) 
    170             CALL vor_een( kt, ncor, ua, va )                ! planetary vorticity trend 
     172            CALL vor_een( kt, ncor, ua, va )                      ! planetary vorticity trend 
    171173            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    172174            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    173175            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 
    174176         ELSE 
    175             CALL vor_een( kt, ntot, ua, va )                ! total vorticity 
     177            CALL vor_een( kt, ntot, ua, va )                ! total vorticity trend 
    176178         ENDIF 
    177179         ! 
     
    197199      !! 
    198200      !! ** Method  :   Trend evaluated using now fields (centered in time)  
    199       !!      and the Sadourny (1975) flux form formulation : conserves the 
    200       !!      horizontal kinetic energy. 
    201       !!      The trend of the vorticity term is given by: 
    202       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    203       !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f  mi(e1v*e3v vn) ] 
    204       !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f  mj(e2u*e3u un) ] 
    205       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    206       !!          voru = 1/e1u  mj-1[ (rotn+f)  mi(e1v vn) ] 
    207       !!          vorv = 1/e2v  mi-1[ (rotn+f)  mj(e2u un) ] 
    208       !!      Add this trend to the general momentum trend (ua,va): 
    209       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
     201      !!       and the Sadourny (1975) flux form formulation : conserves the 
     202      !!       horizontal kinetic energy. 
     203      !!         The general trend of momentum is increased due to the vorticity  
     204      !!       term which is given by: 
     205      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ] 
     206      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ] 
     207      !!       where rvor is the relative vorticity 
    210208      !! 
    211209      !! ** Action : - Update (ua,va) with the now vorticity term trend 
     
    219217      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    220218      ! 
    221       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    222       REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars 
    223       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz 
     219      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     220      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
     221      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz   ! 2D workspace 
    224222      !!---------------------------------------------------------------------- 
    225223      ! 
     
    233231         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    234232      ENDIF 
    235  
    236       zfact2 = 0.5 * 0.5      ! Local constant initialization 
    237  
     233      ! 
    238234      !                                                ! =============== 
    239235      DO jk = 1, jpkm1                                 ! Horizontal slab 
    240236         !                                             ! =============== 
    241237         ! 
    242          ! Potential vorticity and horizontal fluxes 
    243          ! ----------------------------------------- 
    244          SELECT CASE( kvor )      ! vorticity considered 
    245          CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    246          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
    247          CASE ( 3 )                                                ! metric term 
     238         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     239         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     240            zwz(:,:) = ff(:,:)  
     241         CASE ( np_RVO )                           !* relative vorticity 
     242            DO jj = 1, jpjm1 
     243               DO ji = 1, fs_jpim1   ! vector opt. 
     244                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     245                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     246               END DO 
     247            END DO 
     248         CASE ( np_MET )                           !* metric term 
    248249            DO jj = 1, jpjm1 
    249250               DO ji = 1, fs_jpim1   ! vector opt. 
     
    253254               END DO 
    254255            END DO 
    255          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
    256          CASE ( 5 )                                                ! total (coriolis + metric) 
    257             DO jj = 1, jpjm1 
    258                DO ji = 1, fs_jpim1   ! vector opt. 
    259                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    260                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    261                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    262                        &       * 0.5 * r1_e1e2f(ji,jj)                                              & 
    263                        &       ) 
    264                END DO 
    265             END DO 
     256         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     257            DO jj = 1, jpjm1 
     258               DO ji = 1, fs_jpim1   ! vector opt. 
     259                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     260                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     261                     &                   * r1_e1e2f(ji,jj) 
     262               END DO 
     263            END DO 
     264         CASE ( np_CME )                           !* Coriolis + metric 
     265            DO jj = 1, jpjm1 
     266               DO ji = 1, fs_jpim1   ! vector opt. 
     267                  zwz(ji,jj) = ff(ji,jj)                                                                        & 
     268                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     269                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     270                       &     * 0.5 * r1_e1e2f(ji,jj) 
     271               END DO 
     272            END DO 
     273         CASE DEFAULT                                             ! error 
     274            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    266275         END SELECT 
     276         ! 
     277         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     278            DO jj = 1, jpjm1 
     279               DO ji = 1, fs_jpim1   ! vector opt. 
     280                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     281               END DO 
     282            END DO 
     283         ENDIF 
    267284 
    268285         IF( ln_sco ) THEN 
     
    274291            zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    275292         ENDIF 
    276  
    277          ! Compute and add the vorticity term trend 
    278          ! ---------------------------------------- 
     293         !                                   !==  compute and add the vorticity term trend  =! 
    279294         DO jj = 2, jpjm1 
    280295            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    283298               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    284299               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    285                pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    286                pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     300               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     301               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    287302            END DO   
    288303         END DO   
     
    297312 
    298313 
    299    SUBROUTINE vor_mix( kt ) 
    300       !!---------------------------------------------------------------------- 
    301       !!                 ***  ROUTINE vor_mix  *** 
    302       !! 
    303       !! ** Purpose :   Compute the now total vorticity trend and add it to 
    304       !!      the general trend of the momentum equation. 
    305       !! 
    306       !! ** Method  :   Trend evaluated using now fields (centered in time) 
    307       !!      Mixte formulation : conserves the potential enstrophy of a hori- 
    308       !!      zontally non-divergent flow for (rotzu x uh), the relative vor- 
    309       !!      ticity term and the horizontal kinetic energy for (f x uh), the 
    310       !!      coriolis term. the now trend of the vorticity term is given by: 
    311       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 
    312       !!          voru = 1/e1u  mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] 
    313       !!              +1/e1u  mj-1[ f/e3f          mi(e1v*e3v vn) ] 
    314       !!          vorv = 1/e2v  mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ] 
    315       !!              +1/e2v  mi-1[ f/e3f          mj(e2u*e3u un) ] 
    316       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    317       !!          voru = 1/e1u  mj-1(rotn) mj-1[ mi(e1v vn) ] 
    318       !!              +1/e1u  mj-1[ f          mi(e1v vn) ] 
    319       !!          vorv = 1/e2v  mi-1(rotn) mi-1[ mj(e2u un) ] 
    320       !!              +1/e2v  mi-1[ f          mj(e2u un) ] 
    321       !!      Add this now trend to the general momentum trend (ua,va): 
    322       !!          (ua,va) = (ua,va) + ( voru , vorv ) 
    323       !! 
    324       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 
    325       !! 
    326       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    327       !!---------------------------------------------------------------------- 
    328       ! 
    329       INTEGER, INTENT(in) ::   kt   ! ocean timestep index 
    330       ! 
    331       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    332       REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars 
    333       REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      - 
    334       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
    335       !!---------------------------------------------------------------------- 
    336       ! 
    337       IF( nn_timing == 1 )  CALL timing_start('vor_mix') 
    338       ! 
    339       CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )  
    340       ! 
    341       IF( kt == nit000 ) THEN 
    342          IF(lwp) WRITE(numout,*) 
    343          IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' 
    344          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    345       ENDIF 
    346  
    347       zfact1 = 0.5 * 0.25      ! Local constant initialization 
    348       zfact2 = 0.5 * 0.5 
    349  
    350 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) 
    351       !                                                ! =============== 
    352       DO jk = 1, jpkm1                                 ! Horizontal slab 
    353          !                                             ! =============== 
    354          ! 
    355          ! Relative and planetary potential vorticity and horizontal fluxes 
    356          ! ---------------------------------------------------------------- 
    357          IF( ln_sco ) THEN         
    358             IF( ln_dynadv_vec ) THEN 
    359                zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) 
    360             ELSE                        
    361                DO jj = 1, jpjm1 
    362                   DO ji = 1, fs_jpim1   ! vector opt. 
    363                      zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    364                         &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    365                         &       * 0.5 / ( e1e2f (ji,jj) * fse3f(ji,jj,jk) ) 
    366                   END DO 
    367                END DO 
    368             ENDIF 
    369             zwz(:,:) = ff  (:,:)    / fse3f(:,:,jk) 
    370             zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 
    371             zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 
    372          ELSE 
    373             IF( ln_dynadv_vec ) THEN 
    374                zww(:,:) = rotn(:,:,jk) 
    375             ELSE                        
    376                DO jj = 1, jpjm1 
    377                   DO ji = 1, fs_jpim1   ! vector opt. 
    378                      zww(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    379                         &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    380                         &       * 0.5 * r1_e1e2f(ji,jj) 
    381                   END DO 
    382                END DO 
    383             ENDIF 
    384             zwz(:,:) = ff (:,:) 
    385             zwx(:,:) = e2u(:,:) * un(:,:,jk) 
    386             zwy(:,:) = e1v(:,:) * vn(:,:,jk) 
    387          ENDIF 
    388  
    389          ! Compute and add the vorticity term trend 
    390          ! ---------------------------------------- 
    391          DO jj = 2, jpjm1 
    392             DO ji = fs_2, fs_jpim1   ! vector opt. 
    393                zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 
    394                zy2 = ( zwy(ji,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj) 
    395                zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 
    396                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj) 
    397                ! enstrophy conserving formulation for relative vorticity term 
    398                zua = zfact1 * ( zww(ji  ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 ) 
    399                zva =-zfact1 * ( zww(ji-1,jj  ) + zww(ji,jj) ) * ( zx1 + zx2 ) 
    400                ! energy conserving formulation for planetary vorticity term 
    401                zcua = zfact2 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    402                zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    403                ! mixed vorticity trend added to the momentum trends 
    404                ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua 
    405                va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva 
    406             END DO   
    407          END DO   
    408          !                                             ! =============== 
    409       END DO                                           !   End of slab 
    410       !                                                ! =============== 
    411       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )  
    412       ! 
    413       IF( nn_timing == 1 )  CALL timing_stop('vor_mix') 
    414       ! 
    415    END SUBROUTINE vor_mix 
    416  
    417  
    418314   SUBROUTINE vor_ens( kt, kvor, pua, pva ) 
    419315      !!---------------------------------------------------------------------- 
     
    427323      !!      potential enstrophy of a horizontally non-divergent flow. the 
    428324      !!      trend of the vorticity term is given by: 
    429       !!       * s-coordinate (ln_sco=T), the e3. are inside the derivative: 
    430       !!          voru = 1/e1u  mj-1[ (rotn+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
    431       !!          vorv = 1/e2v  mi-1[ (rotn+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    432       !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 
    433       !!          voru = 1/e1u  mj-1[ rotn+f ]  mj-1[ mi(e1v vn) ] 
    434       !!          vorv = 1/e2v  mi-1[ rotn+f ]  mi-1[ mj(e2u un) ] 
     325      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ] 
     326      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ] 
    435327      !!      Add this trend to the general momentum trend (ua,va): 
    436328      !!          (ua,va) = (ua,va) + ( voru , vorv ) 
     
    440332      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    441333      !!---------------------------------------------------------------------- 
    442       ! 
    443334      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    444335      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
     
    447338      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    448339      ! 
    449       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    450       REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars 
    451       REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
     340      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     341      REAL(wp) ::   zuav, zvau   ! local scalars 
     342      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, zwy, zwz, zww   ! 2D workspace 
    452343      !!---------------------------------------------------------------------- 
    453344      ! 
     
    461352         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    462353      ENDIF 
    463  
    464       zfact1 = 0.5 * 0.25      ! Local constant initialization 
    465  
    466 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 
    467354      !                                                ! =============== 
    468355      DO jk = 1, jpkm1                                 ! Horizontal slab 
    469356         !                                             ! =============== 
    470357         ! 
    471          ! Potential vorticity and horizontal fluxes 
    472          ! ----------------------------------------- 
    473          SELECT CASE( kvor )      ! vorticity considered 
    474          CASE ( 1 )   ;   zwz(:,:) =                  ff(:,:)      ! planetary vorticity (Coriolis) 
    475          CASE ( 2 )   ;   zwz(:,:) =   rotn(:,:,jk)                ! relative  vorticity 
    476          CASE ( 3 )                                                ! metric term 
     358         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     359         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     360            zwz(:,:) = ff(:,:)  
     361         CASE ( np_RVO )                           !* relative vorticity 
     362            DO jj = 1, jpjm1 
     363               DO ji = 1, fs_jpim1   ! vector opt. 
     364                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     365                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     366               END DO 
     367            END DO 
     368         CASE ( np_MET )                           !* metric term 
    477369            DO jj = 1, jpjm1 
    478370               DO ji = 1, fs_jpim1   ! vector opt. 
     
    482374               END DO 
    483375            END DO 
    484          CASE ( 4 )   ;   zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) )    ! total (relative + planetary vorticity) 
    485          CASE ( 5 )                                                ! total (coriolis + metric) 
    486             DO jj = 1, jpjm1 
    487                DO ji = 1, fs_jpim1   ! vector opt. 
    488                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    489                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    490                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    491                        &       * 0.5 * r1_e1e2f(ji,jj)                                                & 
    492                        &       ) 
    493                END DO 
    494             END DO 
     376         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     377            DO jj = 1, jpjm1 
     378               DO ji = 1, fs_jpim1   ! vector opt. 
     379                  zwz(ji,jj) = ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     380                     &                      - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     381                     &                   * r1_e1e2f(ji,jj) 
     382               END DO 
     383            END DO 
     384         CASE ( np_CME )                           !* Coriolis + metric 
     385            DO jj = 1, jpjm1 
     386               DO ji = 1, fs_jpim1   ! vector opt. 
     387                  zwz(ji,jj) = ff(ji,jj)                                                                       & 
     388                       &     + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     389                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     390                       &     * 0.5 * r1_e1e2f(ji,jj) 
     391               END DO 
     392            END DO 
     393         CASE DEFAULT                                             ! error 
     394            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    495395         END SELECT 
     396         ! 
     397         IF( ln_dynvor_msk ) THEN           !==  mask/unmask vorticity ==! 
     398            DO jj = 1, jpjm1 
     399               DO ji = 1, fs_jpim1   ! vector opt. 
     400                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     401               END DO 
     402            END DO 
     403         ENDIF 
    496404         ! 
    497405         IF( ln_sco ) THEN                   !==  horizontal fluxes  ==! 
     
    506414         DO jj = 2, jpjm1 
    507415            DO ji = fs_2, fs_jpim1   ! vector opt. 
    508                zuav = zfact1 * r1_e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
    509                   &                            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) 
    510                zvau =-zfact1 * r1_e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    511                   &                            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
     416               zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   & 
     417                  &                       + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) 
     418               zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
     419                  &                       + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
    512420               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    513421               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     
    542450      !!---------------------------------------------------------------------- 
    543451      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    544       INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; 
    545       !                                                           ! =nrvm (relative vorticity or metric) 
     452      INTEGER , INTENT(in   )                         ::   kvor   ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 
    546453      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pua    ! total u-trend 
    547454      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    548       !! 
    549       INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
    550       INTEGER  ::   ierr                                          ! local integer 
    551       REAL(wp) ::   zfac12, zua, zva                              ! local scalars 
    552       REAL(wp) ::   zmsk, ze3                                     ! local scalars 
    553       !                                                           !  3D workspace  
    554       REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz 
    555       REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse 
    556 #if defined key_vvl 
    557       REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T) 
    558 #else 
    559       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all 
    560 #endif 
     455      ! 
     456      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     457      INTEGER  ::   ierr         ! local integer 
     458      REAL(wp) ::   zua, zva     ! local scalars 
     459      REAL(wp) ::   zmsk, ze3    ! local scalars 
     460      ! 
     461      REAL(wp), POINTER, DIMENSION(:,:)   :: zwx, zwy, zwz, z1_e3f 
     462      REAL(wp), POINTER, DIMENSION(:,:)   :: ztnw, ztne, ztsw, ztse 
    561463      !!---------------------------------------------------------------------- 
    562464      ! 
    563465      IF( nn_timing == 1 )  CALL timing_start('vor_een') 
    564466      ! 
    565       CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        )  
    566       CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
    567 #if defined key_vvl 
    568       CALL wrk_alloc( jpi, jpj, jpk, ze3f                   ) 
    569 #endif 
     467      CALL wrk_alloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f )  
     468      CALL wrk_alloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   )  
    570469      ! 
    571470      IF( kt == nit000 ) THEN 
     
    573472         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
    574473         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    575 #if ! defined key_vvl 
    576          IF( .NOT.ALLOCATED(ze3f) ) THEN 
    577             ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr ) 
    578             IF( lk_mpp    )   CALL mpp_sum ( ierr ) 
    579             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' ) 
    580          ENDIF 
    581          ze3f(:,:,:) = 0._wp 
    582 #endif 
    583474      ENDIF 
    584  
    585       IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 
    586  
    587          IF( ln_dynvor_een_old ) THEN ! original formulation 
    588             DO jk = 1, jpk 
    589                DO jj = 1, jpjm1 
    590                   DO ji = 1, jpim1 
    591                      ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    592                         &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    593                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3 
    594                   END DO 
    595                END DO 
    596             END DO 
    597          ELSE ! new formulation from NEMO 3.6 
    598             DO jk = 1, jpk 
    599                DO jj = 1, jpjm1 
    600                   DO ji = 1, jpim1 
    601                      ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    602                         &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    603                      zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    604                         &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
    605                      IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
    606                   END DO 
    607                END DO 
    608             END DO 
    609          ENDIF 
    610  
    611          CALL lbc_lnk( ze3f, 'F', 1. ) 
    612       ENDIF 
    613  
    614       zfac12 = 1._wp / 12._wp    ! Local constant initialization 
    615  
     475      ! 
    616476      !                                                ! =============== 
    617477      DO jk = 1, jpkm1                                 ! Horizontal slab 
    618478         !                                             ! =============== 
    619           
    620          ! Potential vorticity and horizontal fluxes 
    621          ! ----------------------------------------- 
    622          SELECT CASE( kvor )      ! vorticity considered 
    623          CASE ( 1 )                                                ! planetary vorticity (Coriolis) 
    624             zwz(:,:) = ff(:,:)      * ze3f(:,:,jk) 
    625          CASE ( 2 )                                                ! relative  vorticity 
    626             zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 
    627          CASE ( 3 )                                                ! metric term 
     479         ! 
     480         SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point 
     481         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     482            DO jj = 1, jpjm1 
     483               DO ji = 1, fs_jpim1   ! vector opt. 
     484                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     485                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
     486                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4.0_wp / ze3 
     487                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp 
     488                  ENDIF 
     489               END DO 
     490            END DO 
     491         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     492            DO jj = 1, jpjm1 
     493               DO ji = 1, fs_jpim1   ! vector opt. 
     494                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     495                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
     496                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
     497                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
     498                  IF( ze3 /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = zmsk / ze3 
     499                  ELSE                      ;   z1_e3f(ji,jj) = 0.0_wp 
     500                  ENDIF 
     501               END DO 
     502            END DO 
     503         END SELECT 
     504         ! 
     505         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     506         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     507            DO jj = 1, jpjm1 
     508               DO ji = 1, fs_jpim1   ! vector opt. 
     509                  zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 
     510               END DO 
     511            END DO 
     512         CASE ( np_RVO )                           !* relative vorticity 
     513            DO jj = 1, jpjm1 
     514               DO ji = 1, fs_jpim1   ! vector opt. 
     515                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     516                     &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     517                     &       * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     518               END DO 
     519            END DO 
     520         CASE ( np_MET )                           !* metric term 
    628521            DO jj = 1, jpjm1 
    629522               DO ji = 1, fs_jpim1   ! vector opt. 
    630523                  zwz(ji,jj) = (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    631524                       &         - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    632                        &     * 0.5 * r1_e1e2f(ji,jj) * ze3f(ji,jj,jk) 
    633                END DO 
    634             END DO 
    635             CALL lbc_lnk( zwz, 'F', 1. ) 
    636         CASE ( 4 )                                                ! total (relative + planetary vorticity) 
    637             zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 
    638          CASE ( 5 )                                                ! total (coriolis + metric) 
    639             DO jj = 1, jpjm1 
    640                DO ji = 1, fs_jpim1   ! vector opt. 
    641                   zwz(ji,jj) = ( ff (ji,jj)                                                                       & 
    642                        &       + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
    643                        &           - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
    644                        &       * 0.5 * r1_e1e2f(ji,jj)   ) * ze3f(ji,jj,jk) 
    645                END DO 
    646             END DO 
    647             CALL lbc_lnk( zwz, 'F', 1. ) 
     525                       &     * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 
     526               END DO 
     527            END DO 
     528         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     529            DO jj = 1, jpjm1 
     530               DO ji = 1, fs_jpim1   ! vector opt. 
     531                  zwz(ji,jj) = (  ff(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
     532                     &                         - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) & 
     533                     &                      * r1_e1e2f(ji,jj)    ) * z1_e3f(ji,jj) 
     534               END DO 
     535            END DO 
     536         CASE ( np_CME )                           !* Coriolis + metric 
     537            DO jj = 1, jpjm1 
     538               DO ji = 1, fs_jpim1   ! vector opt. 
     539                  zwz(ji,jj) = (  ff(ji,jj)                                                                        & 
     540                       &        + (   ( vn(ji+1,jj  ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj  ) - e2v(ji,jj) )       & 
     541                       &            - ( un(ji  ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji  ,jj+1) - e1u(ji,jj) )   )   & 
     542                       &        * 0.5 * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj) 
     543               END DO 
     544            END DO 
     545         CASE DEFAULT                                             ! error 
     546            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
    648547         END SELECT 
     548         ! 
     549         CALL lbc_lnk( zwz, 'F', 1. ) 
     550         ! 
     551         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     552            DO jj = 1, jpjm1 
     553               DO ji = 1, fs_jpim1   ! vector opt. 
     554                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     555               END DO 
     556            END DO 
     557         ENDIF 
    649558         ! 
    650559         !                                   !==  horizontal fluxes  ==! 
     
    671580         DO jj = 2, jpjm1 
    672581            DO ji = fs_2, fs_jpim1   ! vector opt. 
    673                zua = + zfac12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    674                   &                           + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    675                zva = - zfac12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    676                   &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     582               zua = + r1_12 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     583                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     584               zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     585                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    677586               pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
    678587               pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     
    682591      END DO                                           !   End of slab 
    683592      !                                                ! =============== 
    684       CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        )  
    685       CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
    686 #if defined key_vvl 
    687       CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   ) 
    688 #endif 
     593      ! 
     594      CALL wrk_dealloc( jpi,jpj,   zwx , zwy , zwz , z1_e3f )  
     595      CALL wrk_dealloc( jpi,jpj,   ztnw, ztne, ztsw, ztse   )  
    689596      ! 
    690597      IF( nn_timing == 1 )  CALL timing_stop('vor_een') 
     
    704611      INTEGER ::   ios             ! Local integer output status for namelist read 
    705612      !! 
    706       NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old 
     613      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 
    707614      !!---------------------------------------------------------------------- 
    708615 
     
    721628         WRITE(numout,*) '~~~~~~~~~~~~' 
    722629         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme' 
    723          WRITE(numout,*) '           energy    conserving scheme                ln_dynvor_ene = ', ln_dynvor_ene 
    724          WRITE(numout,*) '           enstrophy conserving scheme                ln_dynvor_ens = ', ln_dynvor_ens 
    725          WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix 
    726          WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een 
    727          WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 
     630         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene 
     631         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens 
     632         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
     633         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
     634         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     635         WRITE(numout,*) '           masked (=1) or unmasked(=0) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
    728636      ENDIF 
    729637 
     638!!gm  this should be removed when choosing a unique strategy for fmask at the coast 
    730639      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 
    731640      ! at angles with three ocean points and one land point 
     641      IF(lwp) WRITE(numout,*) 
     642      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat 
    732643      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 
    733644         DO jk = 1, jpk 
     
    743654          ! 
    744655      ENDIF 
    745  
    746       ioptio = 0                     ! Control of vorticity scheme options 
    747       IF( ln_dynvor_ene )   ioptio = ioptio + 1 
    748       IF( ln_dynvor_ens )   ioptio = ioptio + 1 
    749       IF( ln_dynvor_mix )   ioptio = ioptio + 1 
    750       IF( ln_dynvor_een )   ioptio = ioptio + 1 
    751       IF( ln_dynvor_een_old )   ioptio = ioptio + 1 
    752       IF( lk_esopa      )   ioptio =          1 
    753  
     656!!gm end 
     657 
     658      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme) 
     659      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENE   ;   ENDIF 
     660      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_ENS   ;   ENDIF 
     661      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_MIX   ;   ENDIF 
     662      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF 
     663      ! 
    754664      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
    755  
    756       !                              ! Set nvor (type of scheme for vorticity) 
    757       IF( ln_dynvor_ene )   nvor =  0 
    758       IF( ln_dynvor_ens )   nvor =  1 
    759       IF( ln_dynvor_mix )   nvor =  2 
    760       IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3 
    761       IF( lk_esopa      )   nvor = -1 
    762        
    763       !                              ! Set ncor, nrvm, ntot (type of vorticity) 
    764       IF(lwp) WRITE(numout,*) 
    765       ncor = 1 
     665      !                       
     666      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot) 
     667      ncor = np_COR 
    766668      IF( ln_dynadv_vec ) THEN      
    767669         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity' 
    768          nrvm = 2 
    769          ntot = 4 
     670         nrvm = np_RVO        ! relative vorticity 
     671         ntot = np_CRV        ! relative + planetary vorticity 
    770672      ELSE                         
    771673         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term' 
    772          nrvm = 3 
    773          ntot = 5 
     674         nrvm = np_MET        ! metric term 
     675         ntot = np_CME        ! Coriolis + metric term 
    774676      ENDIF 
    775677       
    776678      IF(lwp) THEN                   ! Print the choice 
    777679         WRITE(numout,*) 
    778          IF( nvor ==  0 )   WRITE(numout,*) '         vorticity scheme : energy conserving scheme' 
    779          IF( nvor ==  1 )   WRITE(numout,*) '         vorticity scheme : enstrophy conserving scheme' 
    780          IF( nvor ==  2 )   WRITE(numout,*) '         vorticity scheme : mixed enstrophy/energy conserving scheme' 
    781          IF( nvor ==  3 )   WRITE(numout,*) '         vorticity scheme : energy and enstrophy conserving scheme' 
    782          IF( nvor == -1 )   WRITE(numout,*) '         esopa test: use all lateral physics options' 
     680         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme' 
     681         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme' 
     682         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 
     683         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme' 
    783684      ENDIF 
    784685      ! 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5737 r5777  
    4949      !! 
    5050      !! ** Method  :   The now vertical advection of momentum is given by: 
    51       !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 
    52       !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 
     51      !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ] 
     52      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ] 
    5353      !!      Add this trend to the general trend (ua,va): 
    5454      !!         (ua,va) = (ua,va) + w dz(u,v) 
     
    183183      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts') 
    184184      ! 
    185       CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww )  
    186       CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs )  
     185      CALL wrk_alloc( jpi,jpj,jpk,     zwuw, zwvw, zww )  
     186      CALL wrk_alloc( jpi,jpj,jpk,3,   zus , zvs )  
    187187      ! 
    188188      IF( kt == nit000 ) THEN 
     
    210210         END DO 
    211211      END DO 
    212       ! 
    213       ! Surface and bottom advective fluxes set to zero 
    214       DO jj = 2, jpjm1         
     212 
     213      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero 
    215214         DO ji = fs_2, fs_jpim1           ! vector opt. 
     215 !!gm missing ISF boundary condition 
    216216            zwuw(ji,jj, 1 ) = 0._wp 
    217217            zwvw(ji,jj, 1 ) = 0._wp 
     
    284284         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    285285      ! 
    286       CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww )  
    287       CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs )  
     286      CALL wrk_dealloc( jpi,jpj,jpk,     zwuw, zwvw, zww )  
     287      CALL wrk_dealloc( jpi,jpj,jpk,3,   zus , zvs )  
    288288      ! 
    289289      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts') 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r5760 r5777  
    99 
    1010   !!---------------------------------------------------------------------- 
    11    !!   dyn_zdf      : Update the momentum trend with the vertical diffusion 
    12    !!   dyn_zdf_init : initializations of the vertical diffusion scheme 
     11   !!   dyn_zdf       : Update the momentum trend with the vertical diffusion 
     12   !!   dyn_zdf_init  : initializations of the vertical diffusion scheme 
    1313   !!---------------------------------------------------------------------- 
    14    USE oce             ! ocean dynamics and tracers variables 
    15    USE dom_oce         ! ocean space and time domain variables  
    16    USE zdf_oce         ! ocean vertical physics variables 
    17  
    18    USE dynzdf_exp      ! vertical diffusion: explicit (dyn_zdf_exp     routine) 
    19    USE dynzdf_imp      ! vertical diffusion: implicit (dyn_zdf_imp     routine) 
    20  
    21    USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    22    USE trd_oce         ! trends: ocean variables 
    23    USE trddyn          ! trend manager: dynamics 
    24    USE in_out_manager  ! I/O manager 
    25    USE lib_mpp         ! MPP library 
    26    USE prtctl          ! Print control 
    27    USE wrk_nemo        ! Memory Allocation 
    28    USE timing          ! Timing 
     14   USE oce            ! ocean dynamics and tracers variables 
     15   USE dom_oce        ! ocean space and time domain variables  
     16   USE zdf_oce        ! ocean vertical physics variables 
     17   USE dynzdf_exp     ! vertical diffusion: explicit (dyn_zdf_exp     routine) 
     18   USE dynzdf_imp     ! vertical diffusion: implicit (dyn_zdf_imp     routine) 
     19   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
     20   USE trd_oce        ! trends: ocean variables 
     21   USE trddyn         ! trend manager: dynamics 
     22   ! 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_mpp        ! MPP library 
     25   USE prtctl         ! Print control 
     26   USE wrk_nemo       ! Memory Allocation 
     27   USE timing         ! Timing 
    2928 
    3029   IMPLICIT NONE 
     
    6160      !!--------------------------------------------------------------------- 
    6261      ! 
    63       IF( nn_timing == 1 )  CALL timing_start('dyn_zdf') 
     62      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf') 
    6463      ! 
    6564      !                                          ! set time step 
     
    7978      CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme 
    8079      ! 
    81       CASE ( -1 )                                        ! esopa: test all possibility with control print 
    82                        CALL dyn_zdf_exp( kt, r2dt ) 
    83                        CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               & 
    84                           &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    85                        CALL dyn_zdf_imp( kt, r2dt ) 
    86                        CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               & 
    87                           &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    8880      END SELECT 
    8981 
     
    9688      !                                          ! print mean trends (used for debugging) 
    9789      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               & 
    98             &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    99       ! 
    100       IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf') 
     90         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     91         ! 
     92      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf') 
    10193      ! 
    10294   END SUBROUTINE dyn_zdf 
     
    126118      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    127119      ! 
    128       IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used 
    129       ! 
    130120      IF(lwp) THEN                                  ! Print the choice 
    131121         WRITE(numout,*) 
    132122         WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 
    133123         WRITE(numout,*) '~~~~~~~~~~~' 
    134          IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used' 
    135124         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme' 
    136125         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme' 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r5120 r5777  
    209209      !----------------------------------------------------------------------- 
    210210      ! 
    211       !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    212       DO jk = 2, jpkm1 
     211      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    213212         DO jj = 2, jpjm1    
    214213            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    309308      !----------------------------------------------------------------------- 
    310309      ! 
    311       !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    312       DO jk = 2, jpkm1         
     310      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    313311         DO jj = 2, jpjm1    
    314312            DO ji = fs_2, fs_jpim1   ! vector opt. 
  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r5770 r5777  
    2020   USE sbc_oce         ! surface boundary condition: ocean 
    2121   USE domvvl          ! Variable volume 
    22    USE divcur          ! hor. divergence and curl      (div & cur routines) 
    23    USE restart         ! only for lrst_oce 
    24    USE in_out_manager  ! I/O manager 
    25    USE prtctl          ! Print control 
    26    USE phycst 
    27    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    28    USE lib_mpp         ! MPP library 
     22   USE divhor          ! horizontal divergence 
     23   USE phycst          ! physical constants 
    2924   USE bdy_oce 
    3025   USE bdy_par          
     
    3631   USE asminc          ! Assimilation increment 
    3732#endif 
     33   USE in_out_manager  ! I/O manager 
     34   USE restart         ! only for lrst_oce 
     35   USE prtctl          ! Print control 
     36   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
     37   USE lib_mpp         ! MPP library 
    3838   USE wrk_nemo        ! Memory Allocation 
    3939   USE timing          ! Timing 
     
    6666      !!      by the time step. 
    6767      !! 
    68       !! ** action  :   ssha    : after sea surface height 
     68      !! ** action  :   ssha, after sea surface height 
    6969      !! 
    7070      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7171      !!---------------------------------------------------------------------- 
    72       ! 
    73       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv 
    74       INTEGER, INTENT(in) ::   kt                      ! time step 
     72      INTEGER, INTENT(in) ::   kt   ! time step 
    7573      !  
    76       INTEGER             ::   jk                      ! dummy loop indice 
    77       REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
    78       !!---------------------------------------------------------------------- 
    79       ! 
    80       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     74      INTEGER  ::   jk            ! dummy loop indice 
     75      REAL(wp) ::   z2dt, zcoef   ! local scalars 
     76      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace 
     77      !!---------------------------------------------------------------------- 
     78      ! 
     79      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt') 
    8180      ! 
    8281      CALL wrk_alloc( jpi,jpj,   zhdiv )  
     
    8887      ENDIF 
    8988      ! 
    90       CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
     89      CALL div_hor( kt )                              ! Horizontal divergence 
    9190      ! 
    9291      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
     
    104103      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    105104      !  
    106       z1_rau0 = 0.5_wp * r1_rau0 
    107       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     105      zcoef = 0.5_wp * r1_rau0 
     106      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    108107 
    109108#if ! defined key_dynspg_ts 
    110109      ! These lines are not necessary with time splitting since 
    111110      ! boundary condition on sea level is set during ts loop 
    112 #if defined key_agrif 
     111# if defined key_agrif 
    113112      CALL agrif_ssh( kt ) 
    114 #endif 
    115 #if defined key_bdy 
    116       IF (lk_bdy) THEN 
    117          CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 
    118          CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 
    119       ENDIF 
    120 #endif 
     113# endif 
     114# if defined key_bdy 
     115      IF( lk_bdy ) THEN <