Changeset 2887


Ignore:
Timestamp:
2011-10-05T11:16:29+02:00 (9 years ago)
Author:
hliu
Message:

addition and modification of files for Smagorinsky method. for Maria Luneva

Location:
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM
Files:
2 added
14 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/CONFIG/GYRE/cpp_GYRE.fcm

    r2670 r2887  
    1  bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke key_vectopt_loop key_iomput 
     1 bld::tool::fppkeys key_gyre key_dynspg_flt key_ldfslp key_zdftke   key_traldf_c3d key_traldf_smag key_dynldf_c3d key_dynldf_smag key_mpp_mpi 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r2715 r2887  
    44   !! Ocean dynamics:  lateral viscosity trend 
    55   !!====================================================================== 
    6    !! History :  OPA  ! 1990-09  (G. Madec)  Original code 
    7    !!            4.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon) 
    8    !!            6.0  ! 1996-01  (G. Madec)  statement function for e3 
    9    !!            8.0  ! 1997-07  (G. Madec)  lbc calls 
    10    !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module 
    11    !!            2.0  ! 2004-08  (C. Talandier) New trends organization 
    12    !!---------------------------------------------------------------------- 
    136 
    147   !!---------------------------------------------------------------------- 
     
    169   !!                   using an iso-level bilaplacian operator 
    1710   !!---------------------------------------------------------------------- 
     11   !! * Modules used 
    1812   USE oce             ! ocean dynamics and tracers 
    1913   USE dom_oce         ! ocean space and time domain 
     
    2721   PRIVATE 
    2822 
    29    PUBLIC   dyn_ldf_bilap   ! called by step.F90 
     23   !! * Routine accessibility 
     24   PUBLIC dyn_ldf_bilap  ! called by step.F90 
    3025 
    3126   !! * Substitutions 
     
    3631   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3732   !! $Id$  
    38    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    39    !!---------------------------------------------------------------------- 
     33   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     34   !!---------------------------------------------------------------------- 
     35 
    4036CONTAINS 
    4137 
     
    7369      !! ** Action : - Update (ua,va) with the before iso-level biharmonic 
    7470      !!               mixing trend. 
    75       !!---------------------------------------------------------------------- 
    76       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    77       USE wrk_nemo, ONLY:   zcu => wrk_2d_1 , zcv => wrk_2d_2   ! 3D workspace 
    78       USE wrk_nemo, ONLY:   zuf => wrk_3d_3 , zut => wrk_3d_4   ! 3D workspace 
    79       USE wrk_nemo, ONLY:   zlu => wrk_3d_5 , zlv => wrk_3d_6 
    80       ! 
    81       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    82       ! 
    83       INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    84       REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar 
    85       !!---------------------------------------------------------------------- 
    86  
    87       IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 3,4,5,6) ) THEN 
    88          CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable')   ;   RETURN 
    89       ENDIF 
    90  
    91       IF( kt == nit000 .AND. lwp ) THEN 
    92          WRITE(numout,*) 
    93          WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator' 
    94          WRITE(numout,*) '~~~~~~~~~~~~~' 
     71      !! 
     72      !! History : 
     73      !!        !  90-09  (G. Madec)  Original code 
     74      !!        !  91-11  (G. Madec) 
     75      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon) 
     76      !!        !  96-01  (G. Madec)  statement function for e3 
     77      !!        !  97-07  (G. Madec)  lbc calls 
     78      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
     79      !!   9.0  !  04-08  (C. Talandier) New trends organization 
     80      !!---------------------------------------------------------------------- 
     81      !! * Arguments 
     82      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
     83 
     84      !! * Local declarations 
     85      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     86      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v ! temporary scalar 
     87      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     88         zcu, zcv                             ! temporary workspace 
     89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
     90         zuf, zut, zlu, zlv                   ! temporary workspace 
     91      !!---------------------------------------------------------------------- 
     92      !!  OPA 8.5, LODYC-IPSL (2002) 
     93      !!---------------------------------------------------------------------- 
     94 
     95      IF( kt == nit000 ) THEN 
     96         IF(lwp) WRITE(numout,*) 
     97         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator' 
     98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
    9599      ENDIF 
    96100 
     
    100104!!$      zlu(:,:,jpk) = 0.e0 
    101105!!$      zlv(:,:,jpk) = 0.e0 
    102       zuf(:,:,:) = 0._wp 
    103       zut(:,:,:) = 0._wp 
    104       zlu(:,:,:) = 0._wp 
    105       zlv(:,:,:) = 0._wp 
     106      zuf(:,:,:) = 0.e0 
     107      zut(:,:,:) = 0.e0 
     108      zlu(:,:,:) = 0.e0 
     109      zlv(:,:,:) = 0.e0 
    106110 
    107111      !                                                ! =============== 
     
    133137            END DO   
    134138         ENDIF 
    135       END DO 
    136       CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions 
    137  
     139      ENDDO 
     140 
     141      ! Boundary conditions on the laplacian  (zlu,zlv) 
     142      CALL lbc_lnk( zlu, 'U', -1. ) 
     143      CALL lbc_lnk( zlv, 'V', -1. ) 
     144          
    138145          
    139146      DO jk = 1, jpkm1 
     
    143150          
    144151         ! Multiply by the eddy viscosity coef. (at u- and v-points) 
     152#if ! defined key_dynldf_smag 
    145153         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk) 
    146154         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk) 
    147           
     155#endif          
    148156         ! Contravariant "laplacian" 
    149157         zcu(:,:) = e1u(:,:) * zlu(:,:,jk) 
     
    187195         ! Bilaplacian 
    188196         ! ----------- 
    189  
     197#if !defined key_dynldf_smag 
    190198         DO jj = 2, jpjm1 
    191199            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    203211            END DO 
    204212         END DO 
    205  
     213#else 
     214 
     215         DO jj = 2, jpjm1 
     216            DO ji = fs_2, fs_jpim1   ! vector opt. 
     217               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) 
     218               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) 
     219               ! horizontal biharmonic diffusive trends 
     220               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   & 
     221                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj) 
     222 
     223               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   & 
     224                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj) 
     225               ! add it to the general momentum trends 
     226               ua(ji,jj,jk) = ua(ji,jj,jk) + zua*fsahmu(ji,jj,jk) 
     227               va(ji,jj,jk) = va(ji,jj,jk) + zva*fsahmv(ji,jj,jk) 
     228            END DO 
     229         END DO 
     230 
     231#endif 
    206232         !                                             ! =============== 
    207233      END DO                                           !   End of slab 
    208234      !                                                ! =============== 
    209       IF( wrk_not_released(2, 1,2)     .OR.   & 
    210           wrk_not_released(3, 3,4,5,6) )   CALL ctl_stop('dyn_ldf_bilap: failed to release workspace arrays') 
    211       ! 
     235 
    212236   END SUBROUTINE dyn_ldf_bilap 
    213237 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90

    r2715 r2887  
    44   !! Ocean dynamics:  lateral viscosity trend 
    55   !!====================================================================== 
    6    !! History :  OPA  !  1997-07  (G. Madec)  Original code 
    7    !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module 
    8    !!            2.0  !  2004-08  (C. Talandier) New trends organization 
    9    !!---------------------------------------------------------------------- 
    106#if defined key_ldfslp   ||   defined key_esopa 
    117   !!---------------------------------------------------------------------- 
     
    1612   !!   ldfguv         :  
    1713   !!---------------------------------------------------------------------- 
     14   !! * Modules used 
    1815   USE oce             ! ocean dynamics and tracers 
    1916   USE dom_oce         ! ocean space and time domain 
    2017   USE ldfdyn_oce      ! ocean dynamics lateral physics 
    2118   USE zdf_oce         ! ocean vertical physics 
     19   USE in_out_manager  ! I/O manager 
    2220   USE trdmod          ! ocean dynamics trends  
    2321   USE trdmod_oce      ! ocean variables trends 
    2422   USE ldfslp          ! iso-neutral slopes available 
    25    USE in_out_manager  ! I/O manager 
    26    USE lib_mpp         ! MPP library 
    2723   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2824   USE prtctl          ! Print control 
     
    3127   PRIVATE 
    3228 
    33    PUBLIC   dyn_ldf_bilapg       ! called by step.F90 
    34  
    35    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv) 
    36    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv) 
     29   !! * Routine accessibility 
     30   PUBLIC dyn_ldf_bilapg ! called by step.F90 
    3731 
    3832   !! * Substitutions 
     
    4236   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4337   !! $Id$  
    44    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    45    !!---------------------------------------------------------------------- 
     38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     39   !!---------------------------------------------------------------------- 
     40 
    4641CONTAINS 
    47  
    48    INTEGER FUNCTION dyn_ldf_bilapg_alloc() 
    49       !!---------------------------------------------------------------------- 
    50       !!               ***  ROUTINE dyn_ldf_bilapg_alloc  *** 
    51       !!---------------------------------------------------------------------- 
    52       ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     & 
    53          &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc ) 
    54          ! 
    55       IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays') 
    56    END FUNCTION dyn_ldf_bilapg_alloc 
    57  
    5842 
    5943   SUBROUTINE dyn_ldf_bilapg( kt ) 
     
    8367      !!                biharmonic mixing trend. 
    8468      !!              - save the trend in (zwk3,zwk4) ('key_trddyn') 
     69      !! 
     70      !! History : 
     71      !!   8.0  !  97-07  (G. Madec)  Original code 
     72      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
     73      !!   9.0  !  04-08  (C. Talandier) New trends organization 
    8574      !!---------------------------------------------------------------------- 
    86       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    87       USE wrk_nemo, ONLY:   zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4   ! 3D workspace 
    88       USE oce     , ONLY:   zwk3 => ta       , zwk4 => sa         ! ta, sa used as 3D workspace    
    89       ! 
     75      !! * Modules used      
     76      USE oce, ONLY :    zwk3 => ta,   & ! use ta as 3D workspace    
     77                         zwk4 => sa      ! use sa as 3D workspace    
     78 
     79      !! * Arguments 
    9080      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    91       ! 
     81 
     82      !! * Local declarations 
    9283      INTEGER ::   ji, jj, jk                 ! dummy loop indices 
     84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
     85         zwk1, zwk2                ! work array used for rotated biharmonic 
     86         !                         ! operator on tracers and/or momentum 
    9387      !!---------------------------------------------------------------------- 
    94  
    95       IF( wrk_in_use(3, 3,4) ) THEN 
    96          CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN 
    97       ENDIF 
    9888 
    9989      IF( kt == nit000 ) THEN 
     
    10393         zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0 
    10494         zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0 
    105          !                                      ! allocate dyn_ldf_bilapg arrays 
    106          IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') 
    10795      ENDIF 
    10896 
    10997      ! Laplacian of (ub,vb) multiplied by ahm 
    11098      ! --------------------------------------   
    111       CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb) 
    112       !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) 
    113       CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions 
     99      ! rotated harmonic operator applied to (ub,vb) 
     100      !     and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) 
     101 
     102      CALL ldfguv ( ub, vb, zwk1, zwk2, 1 ) 
     103 
     104 
     105      ! Lateral boundary conditions on (zwk1,zwk2) 
     106      CALL lbc_lnk( zwk1, 'U', -1. ) 
     107      CALL lbc_lnk( zwk2, 'V', -1. ) 
     108 
    114109 
    115110      ! Bilaplacian of (ub,vb) 
    116111      ! ----------------------  
    117       CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)  
    118       !                                         ! (output in (zwk3,zwk4) ) 
    119  
    120       ! Update the momentum trends 
     112      ! rotated harmonic operator applied to (zwk1,zwk2) (output in (zwk3,zwk4) ) 
     113 
     114      CALL ldfguv ( zwk1, zwk2, zwk3, zwk4, 2 ) 
     115 
     116 
     117      ! Update the momentum trends           (j-slab :   2, jpj-1) 
    121118      ! -------------------------- 
    122       DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends 
     119      !                                                ! =============== 
     120      DO jj = 2, jpjm1                                 !  Vertical slab 
     121         !                                             ! =============== 
    123122         DO jk = 1, jpkm1 
    124123            DO ji = 2, jpim1 
     124               ! add the diffusive trend to the general momentum trends 
    125125               ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk) 
    126126               va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk) 
    127127            END DO 
    128128         END DO 
    129       END DO 
    130       ! 
    131       IF( wrk_not_released(3, 3,4) )   CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays') 
    132       ! 
     129         !                                             ! =============== 
     130      END DO                                           !   End of slab 
     131      !                                                ! =============== 
     132 
    133133   END SUBROUTINE dyn_ldf_bilapg 
    134134 
     
    174174      !!                          second order vertical derivative term) 
    175175      !!      'key_trddyn' defined: the trend is saved for diagnostics. 
     176      !! 
     177      !! History : 
     178      !!   8.0  !  97-07  (G. Madec)  Original code 
     179      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module 
    176180      !!---------------------------------------------------------------------- 
    177       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    178       USE wrk_nemo, ONLY:   ziut => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt  => wrk_2d_3 
    179       USE wrk_nemo, ONLY:   zivf => wrk_2d_4 , zdku  => wrk_2d_5 , zdk1u => wrk_2d_6 
    180       USE wrk_nemo, ONLY:   zdkv => wrk_2d_7 , zdk1v => wrk_2d_8 
    181       !! 
    182       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity  
    183       !                                                               ! 2nd call: ahm x these fields 
    184       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to 
    185       !                                                               ! pu and pv (all the components except 
    186       !                                                               ! second order vertical derivative term) 
    187       INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call 
    188       ! 
    189       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    190       REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar 
    191       REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      - 
    192       REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      - 
    193       REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
     181      !! * Arguments 
     182      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   & 
     183         pu, pv     ! momentum fields (before u and v for the 1st call, and 
     184      !             ! laplacian of these fields multiplied by ahm for the 2nd 
     185      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   & 
     186         plu, plv   ! partial harmonic operator applied to 
     187      !             ! pu and pv (all the components except 
     188      !             ! second order vertical derivative term) 
     189      INTEGER, INTENT( in ) ::   & 
     190         kahm       ! =1 the laplacian is multiplied by the eddy diffusivity coef. 
     191      !             ! =2 no multiplication 
     192 
     193      !! * Local declarations 
     194      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     195      REAL(wp) ::   & 
     196         zabe1, zabe2, zcof1, zcof2,    &  ! temporary scalars 
     197         zcoef0, zcoef3, zcoef4 
     198      REAL(wp) ::   & 
     199         zbur, zbvr, zmkt, zmkf, zuav, zvav,    & 
     200         zuwslpi, zuwslpj, zvwslpi, zvwslpj 
     201      REAL(wp), DIMENSION(jpi,jpj) ::   & 
     202         ziut, zjuf , zjvt, zivf,       &  ! workspace 
     203         zdku, zdk1u, zdkv, zdk1v 
     204      REAL(wp), DIMENSION(jpi,jpk) ::   & 
     205         zfuw, zfvw, zdiu, zdiv,        &  ! workspace 
     206         zdju, zdj1u, zdjv, zdj1v  
    194207      !!---------------------------------------------------------------------- 
    195208 
    196       IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 
    197          CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable')   ;   RETURN 
    198       END IF 
    199209      !                               ! ********** !   ! =============== 
    200210      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab 
     
    412422         ! II.3 Divergence of vertical fluxes added to the horizontal divergence 
    413423         ! --------------------------------------------------------------------- 
    414  
     424#if defined key_dynldf_smag  
    415425         IF( kahm == 1 ) THEN 
    416426            ! multiply the laplacian by the eddy viscosity coefficient 
     
    448458            STOP 'ldfguv' 
    449459         ENDIF 
     460#else 
     461         IF( kahm == 2 ) THEN 
     462            ! multiply the laplacian by the eddy viscosity coefficient 
     463            DO jk = 1, jpkm1 
     464               DO ji = 2, jpim1 
     465                  ! eddy coef. divided by the volume element 
     466                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) 
     467                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) 
     468                  ! vertical divergence 
     469                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1) 
     470                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1) 
     471                  ! harmonic operator applied to (pu,pv) and multiply by ahm 
     472                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur 
     473                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr 
     474               END DO 
     475            END DO 
     476         ELSEIF( kahm == 1 ) THEN 
     477            ! second call, no multiplication 
     478            DO jk = 1, jpkm1 
     479               DO ji = 2, jpim1 
     480                  ! inverse of the volume element 
     481                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) 
     482                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) 
     483                  ! vertical divergence 
     484                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1) 
     485                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1) 
     486                  ! harmonic operator applied to (pu,pv)  
     487                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur 
     488                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr 
     489               END DO 
     490            END DO 
     491         ELSE 
     492            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm 
     493            IF(lwp)WRITE(numout,*) '         We stop' 
     494            STOP 'ldfguv' 
     495         ENDIF 
     496 
     497#endif 
    450498         !                                             ! =============== 
    451499      END DO                                           !   End of slab 
    452500      !                                                ! =============== 
    453  
    454       IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays') 
    455       ! 
    456501   END SUBROUTINE ldfguv 
    457502 
     
    462507CONTAINS 
    463508   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine 
    464       INTEGER, INTENT(in) :: kt 
    465509      WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt 
    466510   END SUBROUTINE dyn_ldf_bilapg 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90

    r2715 r2887  
    3838   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3939   !! $Id$  
    40    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    41    !!---------------------------------------------------------------------- 
     40   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     41   !!---------------------------------------------------------------------- 
     42 
    4243CONTAINS 
    4344 
     
    6263      !!---------------------------------------------------------------------- 
    6364      INTEGER ::   ioptio         ! ??? 
    64       LOGICAL ::   ll_print = .FALSE.    ! Logical flag for printing viscosity coef. 
     65      LOGICAL :: ll_print = .FALSE.    ! Logical flag for printing viscosity coef. 
    6566      !!  
    6667      NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  & 
    6768         &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   & 
    68          &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp 
     69         &                 rn_ahm_0_lap   , rn_ahmb_0      , rn_ahm_0_blp ,cmsmag1,cmsmag2 
    6970      !!---------------------------------------------------------------------- 
    7071 
     
    8586         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0 
    8687         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0_blp    = ', rn_ahm_0_blp 
     88         WRITE(numout,*) '      smagorinsky coefficient for laplacian       cmsmag1    = ' , cmsmag1 
     89         WRITE(numout,*) '      smagorinsky coefficient for bilaplacian     cmsmag2    = ' , cmsmag2 
    8790      ENDIF 
    8891 
     
    128131      ! Lateral eddy viscosity 
    129132      ! ====================== 
     133#if defined key_dynldf_smag && ! defined key_dynldf_c3d 
     134      IF(lwp) WRITE(numout,*) '   key_dynldf_smag can not be used without key_dynldf_c3d' 
     135#endif 
     136 
    130137#if defined key_dynldf_c3d 
    131138      CALL ldf_dyn_c3d( ll_print )   ! ahm = 3D coef. = F( longitude, latitude, depth ) 
     
    206213      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection 
    207214      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1) 
    208       REAL(wp), INTENT(in   ), DIMENSION          (:) ::   pdep       ! depth of the gridpoint (T, U, V, F) 
    209       REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile 
     215      REAL(wp), INTENT(in   ), DIMENSION        (jpk) ::   pdep       ! depth of the gridpoint (T, U, V, F) 
     216      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pah        ! adimensional vertical profile 
    210217      !! 
    211218      INTEGER  ::   jk           ! dummy loop indices 
     
    248255      REAL(wp), INTENT(in   )                         ::   pwam       ! width of inflection 
    249256      REAL(wp), INTENT(in   )                         ::   pbot       ! bottom value (0<pbot<= 1) 
    250       REAL(wp), INTENT(in   ), DIMENSION      (:,:,:) ::   pdep       ! dep of the gridpoint (T, U, V, F) 
    251       REAL(wp), INTENT(inout), DIMENSION      (:,:,:) ::   pah        ! adimensional vertical profile 
     257      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pdep       ! dep of the gridpoint (T, U, V, F) 
     258      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pah        ! adimensional vertical profile 
    252259      !! 
    253260      INTEGER  ::   jk           ! dummy loop indices 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90

    r2715 r2887  
    2323   REAL(wp), PUBLIC ::   rn_ahm_0_blp    =     0._wp   !: lateral bilaplacian eddy viscosity (m4/s) 
    2424   REAL(wp), PUBLIC ::   ahm0, ahmb0, ahm0_blp         !: OLD namelist names 
    25  
     25   REAL(wp), PUBLIC ::   cmsmag1          = 3._wp       !: constant in laplacian Smagorinsky viscosity 
     26   REAL(wp), PUBLIC ::   cmsmag2          = 3._wp       !: constant in bilaplacian Smagorinsky viscosity 
    2627   !                                                                                  !!! eddy coeff. at U-,V-,W-pts [m2/s] 
    2728#if defined key_dynldf_c3d 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2772 r2887  
    3131   USE in_out_manager  ! I/O manager 
    3232   USE prtctl          ! Print control 
     33   USE lib_mpp 
    3334 
    3435   IMPLICIT NONE 
     
    691692            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    692693            !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    693             zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  ) 
    694             zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  ) 
     694            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik)* ABS( zau )  ) 
     695            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav )  ) 
    695696            !                          !- Slope at u- & v-points (uslpml, vslpml) 
    696             uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 
    697             vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 
     697            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik) 
     698            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik) 
    698699            ! 
    699700            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==! 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90

    r2715 r2887  
    6868         &                 ln_traldf_grif , ln_traldf_gdia,                   & 
    6969         &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       & 
    70          &                 rn_slpmax 
     70         &                 rn_slpmax , chsmag 
    7171      !!---------------------------------------------------------------------- 
    7272 
     
    157157      ENDIF 
    158158#endif 
     159#if defined key_traldf_smag && ! defined key_traldf_c3d 
     160       CALL ctl_stop( 'key_traldf_smag can only be used with key_traldf_c3d' ) 
     161#endif 
     162 
    159163      ! 
    160164   END SUBROUTINE ldf_tra_init 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90

    r2715 r2887  
    8989      ENDIF 
    9090 
    91 # if defined key_traldf_eiv 
     91# if defined key_traldf_eiv  
    9292      aeiu(:,:,:) = aeiv0            ! set aeiu = aeiv at u- and v-points, 
    9393      aeiv(:,:,:) = aeiv0            ! and aeiw at w-point 
     
    108108      CALL lbc_lnk( aeiv, 'V', 1. ) 
    109109      CALL lbc_lnk( aeiw, 'W', 1. ) 
    110 # endif 
    111110 
    112111      IF(lwp .AND. ld_print ) THEN 
     
    121120         CALL prihre(aeiw(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    122121      ENDIF 
     122# endif 
    123123      ! 
    124124   END SUBROUTINE ldf_tra_c3d 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90

    r2715 r2887  
    3030   REAL(wp), PUBLIC ::   rn_aeiv_0       = 2000._wp  !: eddy induced velocity coefficient (m2/s) 
    3131   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp   !: slope limit 
    32  
     32   REAL(wp), PUBLIC ::   chsmag          = 1._wp     !: constant in Smagorinsky diffusivity 
    3333   REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names 
    3434   LOGICAL , PUBLIC ::   l_triad_iso     = .FALSE.   !: calculate triads twice 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r2715 r2887  
    2828   USE diaptr          ! poleward transport diagnostics 
    2929   USE trc_oce         ! share passive tracers/Ocean variables 
    30    USE lib_mpp         ! MPP library 
    3130 
    3231   IMPLICIT NONE 
     
    7473      !!               biharmonic mixing trend. 
    7574      !!---------------------------------------------------------------------- 
    76       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    77       USE oce     , ONLY:   ztu  => ua       , ztv  => va                           ! (ua,va) used as workspace 
    78       USE wrk_nemo, ONLY:   zeeu => wrk_2d_1 , zeev => wrk_2d_2 , zlt => wrk_2d_3   ! 2D workspace 
     75      USE oce         , ztu => ua   ! use ua as workspace 
     76      USE oce         , ztv => va   ! use va as workspace 
    7977      !! 
    8078      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
     
    8785      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    8886      REAL(wp) ::  zbtr, ztra       ! local scalars 
     87      REAL(wp), DIMENSION(jpi,jpj) ::   zeeu, zeev, zlt   ! 2D workspace 
    8988      !!---------------------------------------------------------------------- 
    90  
    91       IF( wrk_in_use(2, 1,2,3) ) THEN 
    92          CALL ctl_stop('tra_ldf_bilap: requested workspace arrays unavailable')   ;   RETURN 
    93       ENDIF 
    9489 
    9590      IF( kt == nit000 )  THEN 
     
    128123               END DO 
    129124            ENDIF 
     125#if ! defined key_traldf_smag 
    130126            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
    131127               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    135131               END DO 
    136132            END DO 
     133#else 
     134            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
     135               DO ji = fs_2, fs_jpim1   ! vector opt. 
     136                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     137                  zlt(ji,jj) = zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
     138                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   ) 
     139               END DO 
     140            END DO 
     141#endif 
    137142            CALL lbc_lnk( zlt, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn) 
    138143 
     
    145150               END DO 
    146151            END DO 
     152#if ! defined key_traldf_smag 
    147153            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend 
    148154               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    154160               END DO 
    155161            END DO 
     162#else 
     163            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend 
     164               DO ji = fs_2, fs_jpim1   ! vector opt. 
     165                  ! horizontal diffusive trends 
     166                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     167                  ztra = zbtr * fsahtt(ji,jj,jk)* (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) 
     168                  ! add it to the general tracer trends 
     169                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     170               END DO 
     171            END DO 
     172 
     173#endif 
    156174            !                                              
    157175         END DO                                           ! Horizontal slab 
     
    165183      END DO                                              ! tracer loop 
    166184      !                                                   ! =========== 
    167       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('tra_ldf_bilap: failed to release workspace arrays') 
    168       ! 
    169185   END SUBROUTINE tra_ldf_bilap 
    170186 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90

    r2715 r2887  
    2424   USE diaptr          ! poleward transport diagnostics  
    2525   USE trc_oce         ! share passive tracers/Ocean variables 
    26    USE lib_mpp         ! MPP library 
    2726 
    2827   IMPLICIT NONE 
     
    6665      !!               biharmonic mixing trend. 
    6766      !!---------------------------------------------------------------------- 
    68       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    69       USE wrk_nemo, ONLY:   wk1 => wrk_4d_1 , wk2 => wrk_4d_2     ! 4D workspace 
    70       ! 
    7167      INTEGER         , INTENT(in   )                      ::   kt       ! ocean time-step index 
    7268      CHARACTER(len=3), INTENT(in   )                      ::   cdtype   ! =TRA or TRC (tracer indicator) 
     
    7470      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields 
    7571      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend  
    76       ! 
    77       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
    78       !!---------------------------------------------------------------------- 
    79  
    80       IF( wrk_in_use(4, 1,2) ) THEN 
    81          CALL ctl_stop('tra_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN 
    82       ENDIF 
     72      !! 
     73      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices 
     74      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) ::   wk1, wk2   ! 4D workspace 
     75      !!---------------------------------------------------------------------- 
    8376 
    8477      IF( kt == nit000 )  THEN 
     
    114107         END DO 
    115108      END DO 
    116       ! 
    117       IF( wrk_not_released(4, 1,2) )   CALL ctl_stop('tra_ldf_bilapg : failed to release workspace arrays.') 
    118109      ! 
    119110   END SUBROUTINE tra_ldf_bilapg 
     
    158149      !! 
    159150      !!---------------------------------------------------------------------- 
    160       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz 
    161       USE oce     , ONLY:   zftv => ua       ! ua used as workspace 
    162       USE wrk_nemo, ONLY:   zftu => wrk_2d_1 , zdkt  => wrk_2d_2 , zdk1t => wrk_2d_3 
    163       USE wrk_nemo, ONLY:   zftw => wrk_xz_1 , zdit  => wrk_xz_2  
    164       USE wrk_nemo, ONLY:   zdjt => wrk_xz_3 , zdj1t => wrk_xz_4 
    165       ! 
     151      USE oce         , zftv => ua     ! use ua as workspace 
     152      !! 
    166153      INTEGER         , INTENT(in )                              ::  kt      ! ocean time-step index 
    167154      CHARACTER(len=3), INTENT(in )                              ::  cdtype  ! =TRA or TRC (tracer indicator)  
     
    179166      REAL(wp) ::  zbtr, ztah, ztav 
    180167      REAL(wp) ::  zcof0, zcof1, zcof2, zcof3, zcof4 
    181       !!---------------------------------------------------------------------- 
    182  
    183       IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use_xz(1,2,3,4) )THEN 
    184          CALL ctl_stop('ldfght : requested workspace arrays unavailable')   ;   RETURN 
    185       ENDIF 
     168      REAL(wp), DIMENSION(jpi,jpj) ::  zftu,  zdkt, zdk1t       ! workspace 
     169      REAL(wp), DIMENSION(jpi,jpk) ::  zftw, zdit, zdjt, zdj1t  !  
     170      !!---------------------------------------------------------------------- 
     171 
    186172      ! 
    187173      DO jn = 1, kjpt 
     
    299285 
    300286            ! II.3 Divergence of vertical fluxes added to the horizontal divergence 
    301             ! --------------------------------------------------------------------- 
    302              
     287            ! -------------------------------------------------------------------- 
     288#if ! defined key_traldf_smag             
    303289            IF( kaht == 1 ) THEN 
    304290               ! multiply the laplacian by the eddy diffusivity coefficient 
     
    334320         !                                                ! =============== 
    335321      END DO 
    336       ! 
    337       IF( wrk_not_released(2, 1,2,3)   .OR.   & 
    338           wrk_not_released_xz(1,2,3,4) )   CALL ctl_stop('ldfght : failed to release workspace arrays.') 
    339       ! 
     322#else 
     323            IF( kaht == 2 ) THEN 
     324               ! multiply the laplacian by the eddy diffusivity coefficient 
     325               DO jk = 1, jpkm1 
     326                  DO ji = 2, jpim1   
     327                     ! eddy coef. divided by the volume element 
     328                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     329                     ! vertical divergence 
     330                     ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) ) 
     331                     ! harmonic operator applied to (pt,ps) and multiply by aht 
     332                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 
     333                  END DO 
     334               END DO 
     335            ELSEIF( kaht == 1 ) THEN 
     336               ! second call, no multiplication 
     337               DO jk = 1, jpkm1 
     338                  DO ji = 2, jpim1 
     339                     ! inverse of the volume element 
     340                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     341                     ! vertical divergence 
     342                     ztav = zftw(ji,jk) - zftw(ji,jk+1) 
     343                     ! harmonic operator applied to (pt,ps)  
     344                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 
     345                  END DO 
     346               END DO 
     347            ELSE 
     348               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht 
     349               IF(lwp) WRITE(numout,*) '         We stop' 
     350               STOP 'ldfght' 
     351            ENDIF 
     352            !                                             ! =============== 
     353         END DO                                           !   End of slab 
     354         !                                                ! =============== 
     355      END DO 
     356 
     357 
     358#endif       
    340359   END SUBROUTINE ldfght 
    341360 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step.F90

    r2715 r2887  
    156156      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient 
    157157#endif 
     158#if defined key_traldf_c3d && key_traldf_smag 
     159       CALL ldf_tra_smag( kstp )      ! Smagorinsky diffusivity 
     160#  endif 
     161#if defined key_dynldf_c3d && key_dynldf_smag 
     162       CALL ldf_dyn_smag( kstp )      ! Smagorinsky viscosity 
     163#  endif 
     164 
    158165 
    159166      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r2528 r2887  
    6060   USE ldfslp           ! iso-neutral slopes               (ldf_slp routine) 
    6161   USE ldfeiv           ! eddy induced velocity coef.      (ldf_eiv routine) 
     62   USE ldftra_smag      ! smagorinsky diffusivity   .      (ldf_tra_smag routine) 
     63   USE ldfdyn_smag      ! smagorinsky viscosity            (ldf_dyn_smag routine) 
    6264 
    6365   USE zdftmx           ! tide-induced vertical mixing     (zdf_tmx routine) 
  • branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/TOOLS/COMPILE/cfg.txt

    r2413 r2887  
    1 GYRE OPA_SRC 
    21ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    32GYRE_LOBSTER OPA_SRC TOP_SRC 
     
    65ORCA2_LIM3 OPA_SRC LIM_SRC_3 
    76ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
     7GYRE OPA_SRC 
     8GYRE_5 OPA_SRC 
Note: See TracChangeset for help on using the changeset viewer.