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

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

Changeset 1481


Ignore:
Timestamp:
2009-07-03T17:07:08+02:00 (15 years ago)
Author:
ctlod
Message:

ensure restartability of TKE2 in coupled mode + few minor bugs corrections, see ticket: #466

Location:
trunk/NEMO
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OFF_SRC/IOM/in_out_manager.F90

    r1324 r1481  
    4444   !!---------------------------------------------------------------------- 
    4545   INTEGER            ::   nitrst                 !: time step at which restart file should be written 
    46 #if defined key_zdftke2 
    47    INTEGER            ::   nitrst_tke2            !: time step at which restart file should be written 
    48 #endif 
    4946   !!---------------------------------------------------------------------- 
    5047   !!                    output monitoring 
  • trunk/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r1312 r1481  
    4444   !!---------------------------------------------------------------------- 
    4545   INTEGER            ::   nitrst                 !: time step at which restart file should be written 
    46 #if defined key_zdftke2 
    47    INTEGER            ::   nitrst_tke2            !: time step at which restart file should be written 
    48 #endif 
    4946   !!---------------------------------------------------------------------- 
    5047   !!                    output monitoring 
  • trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r1469 r1481  
    232232         ! 
    233233         !                                           ! Diagnostics and outputs  
    234 #if defined key_zdftke2 
    235          IF( kt .NE. nitend+1 ) THEN 
    236             IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   & 
    237                &             CALL lim_dia  
    238                              CALL lim_wri( 1  )              ! Ice outputs  
    239             IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file  
    240                              CALL lim_var_glo2eqv            ! ??? 
    241             ! 
    242             IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash 
    243             ! 
    244          ENDIF 
    245 #else 
    246234         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   & 
    247235            &             CALL lim_dia  
     
    252240         IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash 
    253241         ! 
    254 #endif 
    255242      ENDIF                                          ! End sea-ice time step only 
    256243 
  • trunk/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

    r1470 r1481  
    171171         ! ---------------- ! 
    172172 
    173 #if defined key_zdftke2 
    174          IF ( kt .NE. nitend+1 ) THEN 
    175173                                        CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file 
    176          ENDIF 
    177 #else 
    178                                         CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file 
    179 #endif 
    180174         IF( .NOT. lk_c1d ) THEN                                        ! Ice dynamics & transport (not in 1D case) 
    181175                                        CALL lim_dyn_2      ( kt )           ! Ice dynamics    ( rheology/dynamics ) 
     
    192186                                        CALL lim_thd_2      ( kt )      ! Ice thermodynamics  
    193187                                        CALL lim_sbc_2      ( kt )      ! Ice/Ocean Mass & Heat fluxes  
    194 #if defined key_zdftke2 
    195          IF( kt .NE. nitend+1 ) THEN 
    196             IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   & 
    197                &                        CALL lim_dia_2      ( kt )      ! Ice Diagnostics  
    198                                         CALL lim_wri_2      ( kt )      ! Ice outputs  
    199             IF( lrst_ice )              CALL lim_rst_write_2( kt )      ! Ice restart file  
    200             ! 
    201          ENDIF 
    202 #else 
     188 
    203189         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   & 
    204190            &                           CALL lim_dia_2      ( kt )      ! Ice Diagnostics 
     
    206192         IF( lrst_ice )                 CALL lim_rst_write_2( kt )      ! Ice restart file 
    207193         ! 
    208 #endif 
    209194      ENDIF 
    210195      ! 
  • trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90

    r1268 r1481  
    5050 
    5151   PUBLIC   zdf_tke2    ! routine called in step module 
     52   PUBLIC   tke2_rst    ! routine called in step module 
    5253 
    5354   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke2 = .TRUE.  !: TKE vertical mixing flag 
    5455   REAL(wp), PUBLIC                         ::   eboost              !: multiplicative coeff of the shear product. 
    5556   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avmt                !: now mixing coefficient of diffusion for TKE 
     57   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   avm                 !: now mixing coefficient of diffusion for TKE 
    5758   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   dissl               !: now mixing lenght of dissipation 
    5859 
     
    166167      INTEGER, INTENT(in) ::   kt ! ocean time step 
    167168 
    168       IF( kt == nit000  ) THEN 
    169          CALL zdf_tke2_init        ! Initialization (first time-step only) 
    170          IF( .NOT. ln_rstart )  CALL tke_tke( kt ) 
    171          CALL tke_mlmc 
    172       ELSE 
    173          CALL tke_tke( kt ) 
    174          IF( kt /= nitend+1 )   CALL tke_mlmc 
    175          IF( lrst_oce ) CALL tke2_rst( kt, 'WRITE' )      ! write en in restart file  
    176       ENDIF 
     169      IF( kt == nit000 )   CALL zdf_tke2_init             ! Initialization (first time-step only) 
     170      ! 
     171      CALL tke_tke                                        ! now tke (en) 
     172      CALL tke_mlmc                                       ! now avmu, avmv, avt 
    177173 
    178174   END SUBROUTINE zdf_tke2 
    179175 
    180    SUBROUTINE tke_tke( kt ) 
     176   SUBROUTINE tke_tke 
    181177      !!---------------------------------------------------------------------- 
    182178      !! Now Turbulent kinetic energy (output in en) 
     
    185181      !! Resolution of a tridiagonal linear system by a "methode de chasse" 
    186182      !! computation from level 2 to jpkm1  (e(1) computed and e(jpk)=0 ). 
    187       !! avmt represents the turbulent coefficient of the TKE 
     183      !! avm represents the turbulent coefficient of the TKE 
    188184      !!---------------------------------------------------------------------- 
    189185      USE oce,     zwd     =>   ua   ! use ua as workspace 
     
    191187      USE oce,     zdiag2  =>   ta   ! use ta as workspace 
    192188      USE oce,     ztkelc  =>   sa   ! use sa as workspace 
    193       ! 
    194       INTEGER, INTENT(in) ::   kt ! ocean time step 
    195189      ! 
    196190      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
     
    209203      zfact1 = -.5 * rdt * rn_efave 
    210204      zfact2 = 1.5 * rdt * rn_ediss 
    211       zfact3 = 0.5 * rdt * rn_ediss 
     205      zfact3 = 0.5       * rn_ediss 
     206 
     207      ! Surface boundary condition on tke 
     208      ! ------------------------------------------------- 
     209      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0) 
     210!CDIR NOVERRCHK 
     211      DO jj = 2, jpjm1 
     212!CDIR NOVERRCHK 
     213         DO ji = fs_2, fs_jpim1   ! vector opt. 
     214            ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
     215            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
     216            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
     217            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 
     218         END DO 
     219      END DO 
    212220 
    213221      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    217225         ! 
    218226         ! Computation of total energy produce by LC : cumulative sum over jk 
    219          zpelc(:,:,1) =  MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 
     227         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 
    220228         DO jk = 2, jpk 
    221             zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 
     229            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 
    222230         END DO 
    223231         ! 
     
    246254# endif 
    247255         ! 
    248          ! TKE Langmuir circulation source term 
     256         ! TKE Langmuir circulation source term added to en 
    249257!CDIR NOVERRCHK 
    250258         DO jk = 2, jpkm1 
     
    259267                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    260268                  ! TKE Langmuir circulation source term 
    261                   ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     269                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) ) * tmask(ji,jj,jk) 
    262270               END DO 
    263271            END DO 
     
    266274      ENDIF 
    267275 
    268       ! Surface boundary condition on tke 
    269       ! ------------------------------------------------- 
    270       ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0) 
    271 !CDIR NOVERRCHK 
    272       DO jj = 2, jpjm1 
    273 !CDIR NOVERRCHK 
    274          DO ji = fs_2, fs_jpim1   ! vector opt. 
    275             ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    276             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    277             zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
    278             en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 
    279          END DO 
    280       END DO 
    281  
    282       ! Preliminary computing 
     276 
     277      ! Shear production at uw- and vw-points (energy conserving form) 
    283278      DO jk = 2, jpkm1 
    284279         DO jj = 2, jpjm1 
    285280            DO ji = fs_2, fs_jpim1   ! vector opt. 
    286281               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) * & 
    287                   &                              ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) 
     282                  &                              ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) / &  
     283                  &                              ( fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) ) 
    288284               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) * & 
    289                   &                              ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) 
     285                  &                              ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) / & 
     286                  &                              ( fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) ) 
    290287            ENDDO 
    291288         ENDDO 
     
    304301         DO jj = 2, jpjm1 
    305302            DO ji = fs_2, fs_jpim1   ! vector opt. 
    306                !                                                             ! shear prod. - stratif. destruction 
    307                zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) 
    308                ! shear 
    309                zdku = zcoef * ( avmu(ji-1,jj  ,jk) + avmu(ji,jj,jk) )        ! shear 
    310                zdkv = zcoef * ( avmv(ji  ,jj-1,jk) + avmv(ji,jj,jk) ) 
    311                zesh2 =  eboost * ( zdku + zdkv ) - avt(ji,jj,jk) * rn2(ji,jj,jk)        ! coefficient (zesh2) 
    312                ! 
     303               !                                                             ! shear prod. at w-point weightened by mask 
     304               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     305                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    313306               !                                                             ! Matrix 
    314307               zcof = zfact1 * tmask(ji,jj,jk) 
    315308               !                                                                  ! lower diagonal 
    316                zdiag2(ji,jj,jk) = zcof * ( avmt (ji,jj,jk  ) + avmt (ji,jj,jk-1) )   & 
     309               zdiag2(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm (ji,jj,jk-1) )   & 
    317310                  &                    / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
    318311               !                                                                  ! upper diagonal 
    319                zdiag1(ji,jj,jk) = zcof * ( avmt (ji,jj,jk+1) + avmt (ji,jj,jk  ) )   & 
     312               zdiag1(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm (ji,jj,jk  ) )   & 
    320313                  &                    / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) ) 
    321314               !                                                                  ! diagonal 
     
    323316               ! 
    324317               !                                                             ! right hand side in en 
    325                IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells 
    326                   en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) & 
    327                      &                        +   rdt  * zesh2 * tmask(ji,jj,jk) 
    328                ELSE                                                               ! Langmuir cell source term 
    329                   en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) & 
    330                      &                        +   rdt  * zesh2 * tmask(ji,jj,jk)                           & 
    331                      &                        +   rdt  * ztkelc(ji,jj,jk) * tmask(ji,jj,jk) 
    332                ENDIF 
     318               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
     319                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk) 
    333320            END DO 
    334321         END DO 
     
    434421      USE oce,     zmxlm  =>   va   ! use va as workspace 
    435422      USE oce,     zmxld  =>   ta   ! use ta as workspace 
    436       USE oce,     visco  =>   sa   ! use sa as workspace 
    437423      ! 
    438424      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
    439       REAL(wp) ::   zbbrau, zrn2, zraug             ! temporary scalars 
    440       REAL(wp) ::   zfact1, ztx2, zdku              !    -         - 
    441       REAL(wp) ::   zfact2, zty2, zdkv              !    -         - 
    442       REAL(wp) ::   zfact3, zcoef, zav              !    -         - 
     425      REAL(wp) ::   zrn2, zraug                     ! temporary scalars 
     426      REAL(wp) ::   ztx2, zdku                      !    -         - 
     427      REAL(wp) ::   zty2, zdkv                      !    -         - 
     428      REAL(wp) ::   zcoef, zav                      !    -         - 
    443429      REAL(wp) ::   zpdl, zri, zsqen                !    -         - 
    444430      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         - 
     
    569555      ! II  Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt ! 
    570556      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<! 
    571       ! Surface boundary condition on avt and avmt : jk = 1 
     557      ! Surface boundary condition on avt and avm : jk = 1 
    572558      ! ------------------------------------------------- 
    573       ! avt (1) = avmb(1) and avt (jpk) = 0. 
    574       ! avmt(1) = avmb(1) and avmt(jpk) = 0. 
     559      ! avt(1) = avmb(1) and avt(jpk) = 0. 
     560      ! avm(1) = avmb(1) and avm(jpk) = 0. 
    575561      ! 
    576562!CDIR NOVERRCHK 
     
    582568               zsqen = SQRT( en(ji,jj,jk) ) 
    583569               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    584                avmt (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 
    585                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    586                visco(ji,jj,jk) = avmt (ji,jj,jk)  
     570               avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk) 
     571               avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    587572               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    588573            END DO 
     
    591576 
    592577      ! Lateral boundary conditions (sign unchanged) 
    593       CALL lbc_lnk( avmt, 'W', 1. )   ;   CALL lbc_lnk( visco, 'W', 1. ) 
     578      CALL lbc_lnk( avm, 'W', 1. ) 
    594579 
    595580      DO jk = 2, jpkm1                                 ! only vertical eddy viscosity 
    596581         DO jj = 2, jpjm1 
    597582            DO ji = fs_2, fs_jpim1   ! vector opt. 
    598                avmu(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji+1,jj  ,jk) ) * umask(ji,jj,jk)   & 
    599                &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) ) 
    600                avmv(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   & 
    601                &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) ) 
     583               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
     584               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
    602585            END DO 
    603586         END DO 
    604587      END DO 
    605588      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    606       ! 
    607       ! 
    608       DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity 
    609          avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk) 
    610          avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk) 
    611       END DO 
     589 
    612590 
    613591      ! Prandtl number 
     
    623601                  zdkv = avmv(ji,jj-1,jk) * (vn(ji,jj-1,jk-1)-vn(ji,jj-1,jk)) * (vb(ji,jj-1,jk-1)-vb(ji,jj-1,jk)) + & 
    624602                    &    avmv(ji,jj  ,jk) * (vn(ji,jj  ,jk-1)-vn(ji,jj  ,jk)) * (vb(ji,jj  ,jk-1)-vb(ji,jj  ,jk)) 
    625                   zri = MAX( rn2(ji,jj,jk), 0. ) / (zcoef * (zdku + zdkv + 1.e-20) / avmt(ji,jj,jk)) ! local Richardson number 
     603                  zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk)/ (zcoef * (zdku + zdkv + 1.e-20) ) ! local Richardson number 
    626604# if defined key_cfg_1d 
    627605                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
     
    770748      DO jk = 1, jpk 
    771749         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
     750         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    772751         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    773752         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
    774          avmt(:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    775753      END DO 
    776754      dissl(:,:,:) = 1.e-12 
     
    796774     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    797775     ! 
    798      INTEGER ::   jit   ! dummy loop indices 
     776     INTEGER ::   jit, jk   ! dummy loop indices 
     777     INTEGER ::   id1, id2, id3, id4, id5, id6 
    799778     !!---------------------------------------------------------------------- 
    800779     ! 
    801      IF( TRIM(cdrw) == 'READ' ) THEN 
    802         IF( ln_rstart ) THEN 
    803            IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN  
     780     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     781        !                                   ! --------------- 
     782        IF( ln_rstart ) THEN                   !* Read the restart file 
     783           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. ) 
     784           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. ) 
     785           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. ) 
     786           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 
     787           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 
     788           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 
     789           ! 
     790           IF( id1 > 0 ) THEN                       ! 'en' exists 
    804791              CALL iom_get( numror, jpdom_autoglo, 'en', en ) 
    805            ELSE 
    806               IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 )   & 
    807                  &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme' 
    808               IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file' 
    809               IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop' 
     792              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist 
     793                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   ) 
     794                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   ) 
     795                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  ) 
     796                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  ) 
     797                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 
     798              ELSE                                       ! one at least array is missing 
     799                 CALL tke_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
     800              ENDIF 
     801           ELSE                                     ! No TKE array found: initialisation 
     802              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 
    810803              en (:,:,:) = rn_emin * tmask(:,:,:) 
    811               DO jit = 2, nn_itke + 1 
    812                  CALL zdf_tke2( jit ) 
    813               END DO 
     804              CALL tke_mlmc                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
     805              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke2( jit )   ;   END DO 
    814806           ENDIF 
    815         ELSE 
    816            en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value 
     807        ELSE                                   !* Start from rest 
     808           en(:,:,:) = rn_emin * tmask(:,:,:) 
     809           DO jk = 1, jpk                           ! set the Kz to the background value 
     810              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
     811              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
     812              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
     813              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     814           END DO 
    817815        ENDIF 
    818      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
     816        ! 
     817     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     818        !                                   ! ------------------- 
    819819        IF(lwp) WRITE(numout,*) '---- tke2-rst ----' 
    820         CALL iom_rstput( kt, nitrst_tke2, numrow, 'en', en ) 
     820        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
     821        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt    ) 
     822        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm    ) 
     823        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu   ) 
     824        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv   ) 
     825        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  ) 
     826        ! 
    821827     ENDIF 
    822828     ! 
  • trunk/NEMO/OPA_SRC/restart.F90

    r1473 r1481  
    2323   USE zpshde          ! partial step: hor. derivative (zps_hde routine) 
    2424   USE eosbn2          ! equation of state            (eos bn2 routine) 
     25   USE zdfddm          ! double diffusion mixing  
     26   USE zdfmxl          ! mixed layer depth 
    2527   USE trdmld_oce      ! ocean active mixed layer tracers trends variables 
    26 #if defined key_zdftke2 
    27    USE zdf_oce 
    28 #endif 
    2928 
    3029   IMPLICIT NONE 
     
    6766         lrst_oce = .FALSE.    
    6867         nitrst = nitend 
    69 #if defined key_zdftke2 
    70          nitrst_tke2 = nitrst + 1 
    71 #endif 
    7268      ENDIF 
    7369      IF( MOD( kt - 1, nstock ) == 0 ) THEN    
     
    7672         IF( nitrst > nitend )   nitrst = nitend   ! make sure we write a restart at the end of the run 
    7773      ENDIF 
    78 #if defined key_zdftke2 
    79       IF ( nitrst_tke2 .NE. kt ) nitrst_tke2 = nitrst + 1 
    80 #endif 
    8174      ! to get better performances with NetCDF format: 
    8275      ! we open and define the ocean restart file one time step before writing the data (-> at nitrst - 1) 
     
    119112      !!---------------------------------------------------------------------- 
    120113 
    121 #if defined key_zdftke2 
    122       IF( kt == nitrst_tke2 ) THEN 
    123          CALL iom_close( numrow )     ! close the restart file (only at last time step) 
    124          IF( .NOT. lk_trdmld )   lrst_oce = .FALSE. 
    125       ELSE 
    126 #endif 
    127114      !                                                                     ! the begining of the run [s] 
    128115      CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt               )   ! dynamics time step 
     
    138125      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un      ) 
    139126      CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn      ) 
    140       IF( lk_vvl ) CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn      ) 
    141127      CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tn      ) 
    142128      CALL iom_rstput( kt, nitrst, numrow, 'sn'     , sn      ) 
    143129      CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    ) 
    144130      CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   ) 
    145       CALL iom_rstput( kt, nitrst, numrow, 'rn2' , rn2  ) 
    146  
    147 #if defined key_zdftke2 
    148          CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop    ) 
    149 #endif 
    150       IF( nn_dynhpg_rst == 1 .OR. lk_vvl ) THEN 
    151          CALL iom_rstput( kt, nitrst, numrow, 'rhd' , rhd  ) 
    152 #if defined key_zdftke2 
    153          CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt  ) 
    154       ENDIF 
    155 #else 
    156          CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop ) 
    157 #endif 
    158          IF( ln_zps ) THEN 
    159             CALL iom_rstput( kt, nitrst, numrow, 'gtu' , gtu ) 
    160             CALL iom_rstput( kt, nitrst, numrow, 'gsu' , gsu ) 
    161             CALL iom_rstput( kt, nitrst, numrow, 'gru' , gru ) 
    162             CALL iom_rstput( kt, nitrst, numrow, 'gtv' , gtv ) 
    163             CALL iom_rstput( kt, nitrst, numrow, 'gsv' , gsv ) 
    164             CALL iom_rstput( kt, nitrst, numrow, 'grv' , grv ) 
    165          ENDIF 
    166       ENDIF 
    167 #if ! defined key_zdftke2 
     131 
     132      CALL iom_rstput( kt, nitrst, numrow, 'rhop'  , rhop   ) 
     133 
     134      IF( lk_zdfddm )   CALL iom_rstput( kt, nitrst, numrow, 'rrau' , rrau  ) 
     135 
    168136      IF( kt == nitrst ) THEN 
    169137         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
    170138         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE. 
    171139      ENDIF 
    172 #endif 
    173140      ! 
    174141   END SUBROUTINE rst_write 
     
    232199      CALL iom_get( numror, jpdom_autoglo, 'un'   , un    )        ! now    i-component velocity 
    233200      CALL iom_get( numror, jpdom_autoglo, 'vn'   , vn    )        ! now    j-component velocity 
    234       IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'wn'   , wn    )        ! now    k-component velocity 
    235201      CALL iom_get( numror, jpdom_autoglo, 'tn'   , tn    )        ! now    temperature 
    236202      CALL iom_get( numror, jpdom_autoglo, 'sn'   , sn    )        ! now    salinity 
     
    238204      CALL iom_get( numror, jpdom_autoglo, 'hdivn', hdivn )        ! now    horizontal divergence 
    239205 
     206      CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop )          ! now    mixed layer depth 
     207 
     208      IF( lk_zdfddm ) THEN 
     209         IF( iom_varid( numror, 'rrau', ldstop = .FALSE. ) > 0 ) THEN 
     210            CALL iom_get( numror, jpdom_autoglo, 'rrau' , rrau  ) 
     211         ELSE 
     212            CALL eos_init                   ! read equation state type neos parameter 
     213            CALL eos( tb, sb, rhd, rhop )   ! compute rrau 
     214         ENDIF 
     215      ENDIF  
    240216 
    241217      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0) 
     
    247223         hdivb(:,:,:) = hdivn(:,:,:) 
    248224      ENDIF 
    249       CALL eos_init            ! usefull to get the equation state type neos parameter 
    250       IF( iom_varid( numror, 'rn2', ldstop = .FALSE. ) > 0 ) THEN 
    251          CALL iom_get( numror, jpdom_autoglo, 'rn2' , rn2  ) 
    252       ELSE 
    253          CALL bn2( tb, sb, rn2 )  ! before Brunt-Vaisala frequency    
    254       ENDIF 
    255  
    256 #if defined key_zdftke2 
    257       CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities 
    258       IF( iom_varid( numror, 'rhd', ldstop = .FALSE. ) > 0 ) THEN 
    259          CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  ) 
    260       ENDIF 
    261       IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN 
    262          CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop ) 
    263       ENDIF 
    264 #else 
    265       IF( iom_varid( numror, 'rhd', ldstop = .FALSE. ) > 0 ) THEN 
    266          CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  ) 
    267          CALL iom_get( numror, jpdom_autoglo, 'rhop', rhop ) 
    268       ELSE 
    269          CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities 
    270       ENDIF 
    271 #endif 
    272 #if defined key_zdftke2 
    273       IF( iom_varid( numror, 'avt', ldstop = .FALSE. ) > 0 ) THEN 
    274          CALL iom_get( numror, jpdom_autoglo, 'avt' , avt  ) 
    275       ELSE 
    276          IF ( ln_dynhpg_imp ) avt (:,:,:) = 1.2e-5 * tmask(:,:,:) 
    277       ENDIF 
    278 #endif 
    279  
    280       IF( ln_zps .AND. .NOT. lk_c1d ) THEN 
    281          IF( iom_varid( numror, 'gtu', ldstop = .FALSE. ) > 0 ) THEN 
    282             CALL iom_get( numror, jpdom_autoglo, 'gtu' , gtu ) 
    283             CALL iom_get( numror, jpdom_autoglo, 'gsu' , gsu ) 
    284             CALL iom_get( numror, jpdom_autoglo, 'gru' , gru ) 
    285             CALL iom_get( numror, jpdom_autoglo, 'gtv' , gtv ) 
    286             CALL iom_get( numror, jpdom_autoglo, 'gsv' , gsv ) 
    287             CALL iom_get( numror, jpdom_autoglo, 'grv' , grv ) 
    288          ELSE 
    289             CALL zps_hde( nit000, tb , sb , rhd,   &  ! Partial steps: before Horizontal DErivative 
    290                &                  gtu, gsu, gru,   &  ! of t, s, rd at the bottom ocean level 
    291                &                  gtv, gsv, grv ) 
    292          ENDIF 
    293       ENDIF 
    294225      ! 
    295226   END SUBROUTINE rst_read 
  • trunk/NEMO/OPA_SRC/step.F90

    r1480 r1481  
    207207      ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    208208      !----------------------------------------------------------------------- 
    209       IF( neuler == 0 .AND. kstp == nit000 ) THEN 
    210                         CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency 
    211                         rn2b(:,:,:) = rn2(:,:,:) 
    212       ELSE 
    213                         rn2b(:,:,:) = rn2(:,:,:)             ! before Brunt-Vaisala frequency   
    214                         CALL bn2( tn, sn, rn2 )              ! now    Brunt-Vaisala frequency 
    215       ENDIF 
    216  
    217 #if defined key_zdftke2 
    218       IF ( ln_dynhpg_imp ) THEN 
    219       !----------------------------------------------------------------------- 
    220       !  LATERAL PHYSICS 
    221       !----------------------------------------------------------------------- 
    222                                CALL zdf_mxl( kstp )                 ! mixed layer depth 
    223          IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing 
    224 #  if defined key_traldf_c2d 
    225          IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient 
    226 #  endif 
    227       ENDIF 
    228 #endif 
     209 
     210                       CALL bn2( tb, sb, rn2b )              ! before Brunt-Vaisala frequency 
     211                       CALL bn2( tn, sn, rn2  )              ! now    Brunt-Vaisala frequency 
     212 
    229213      !----------------------------------------------------------------------- 
    230214      !  VERTICAL PHYSICS 
     
    260244                        CALL zdf_mxl( kstp )                 ! mixed layer depth 
    261245 
    262 #if defined key_zdftke2 
    263       IF( .NOT. ln_dynhpg_imp ) THEN 
    264                         CALL eos( tb, sb, rhd, rhop )        ! now (swap=before) in situ density for dynhpg module 
    265 #endif 
     246      IF( lrst_oce .AND. lk_zdftke2 )   &                    ! write tke2 information in the restart file 
     247         &               CALL tke2_rst( kstp, 'WRITE' ) 
     248 
    266249      !----------------------------------------------------------------------- 
    267250      !  LATERAL PHYSICS 
    268251      !----------------------------------------------------------------------- 
    269       IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing 
     252      IF( lk_ldfslp )   THEN  
     253                           CALL eos( tb, sb, rhd, rhop )          ! before in situ density 
     254         IF( ln_zps )      CALL zps_hde( kstp, tb, sb, rhd,  &    ! Partial steps: before horizontal gradient 
     255            &                                gtu, gsu, gru,  &    ! of t, s, rd at the last ocean level 
     256            &                                gtv, gsv, grv ) 
     257                               CALL ldf_slp( kstp, rhd, rn2b )      ! before slope of the lateral mixing 
     258      ENDIF 
    270259#if defined key_traldf_c2d 
    271260      IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient 
    272261#  endif 
    273 #if defined key_zdftke2 
    274       ENDIF 
    275 #endif 
    276262 
    277263#if defined key_top 
     
    307293                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields 
    308294 
    309 #if ! defined key_zdftke2 
    310       IF( ln_zdfnpc      )   CALL tra_npc    ( kstp )       ! update after fields by non-penetrative convection 
    311                              CALL tra_nxt    ( kstp )       ! tracer fields at next time step 
    312  
    313       IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg 
    314                                CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module 
    315          IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient 
    316             &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level 
    317             &                                     gtv, gsv, grv ) 
    318       ELSE                                                  ! centered hpg (default case) 
    319                                CALL eos( tb, sb, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module 
    320          IF( ln_zps    )       CALL zps_hde( kstp, tb, sb, rhd,   &   ! Partial steps: now horizontal gradient 
    321             &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level 
    322             &                                     gtv, gsv, grv ) 
    323       ENDIF 
    324 #else 
    325295      IF( .NOT. ln_dynhpg_imp  ) THEN                       ! centered hpg (default case) 
    326296                               CALL eos( tn, sn, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module 
     
    333303      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg 
    334304                               CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module 
    335          IF( lk_ldfslp     )   CALL bn2( ta, sa, rn2 )                ! Time-filtered Brunt-Vaisala frequency 
    336305         IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient 
    337306            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level 
    338307            &                                     gtv, gsv, grv ) 
    339308      ENDIF 
    340 #endif 
     309 
    341310      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    342311      ! Dynamics 
     
    400369 
    401370      IF( lk_cpl )   CALL sbc_cpl_snd( kstp )                 ! coupled mode : field exchanges 
    402  
    403 #if defined key_zdftke2 
    404       IF( ( kstp == nitend ).AND.( lrst_oce ) )  THEN 
    405  
    406          CALL day( kstp+1 )             ! Calendar 
    407          !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    408          ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    409          !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    410                             CALL sbc    ( kstp+1 )         ! Sea Boundary Condition (including sea-ice) 
    411          !----------------------------------------------------------------------- 
    412          !  VERTICAL PHYSICS 
    413          !----------------------------------------------------------------------- 
    414          ! N.B. ua, va, ta, sa arrays are used as workspace in this section 
    415          !----------------------------------------------------------------------- 
    416                             CALL bn2( tn, sn, rn2 )            ! now Brunt-Vaisala frequency 
    417          !                                                     ! Vertical eddy viscosity and diffusivity coefficients 
    418          IF( lk_zdftke2 )   CALL zdf_tke2 ( kstp+1 )           ! TKE2 closure scheme for Kz 
    419                             CALL rst_write( kstp+1 )           ! close the restart file 
    420       ENDIF 
    421 #endif 
     371      ! 
    422372      ! 
    423373   END SUBROUTINE stp 
Note: See TracChangeset for help on using the changeset viewer.