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 9863 – NEMO

Changeset 9863


Ignore:
Timestamp:
2018-06-30T12:51:02+02:00 (7 years ago)
Author:
gm
Message:

#1911 (ENHANCE-04): simplified implementation of the Euler stepping at nit000

Location:
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF
Files:
36 edited
1 moved

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/cfgs/SHARED/namelist_ref

    r9838 r9863  
    4242   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    4343   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    44       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
    45       nn_rstctl   =    0      !  restart control ==> activated only if ln_rstart=T 
     44      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
     45      nn_rstctl    =    0     !  restart control ==> activated only if ln_rstart=T 
    4646      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    4747      !                          !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_oce_update.F90

    r9780 r9863  
    1212   !!            3.6  !  2014-09  (R. Benshila)  
    1313   !!---------------------------------------------------------------------- 
     14 
     15   !!---------------------------------------------------------------------- 
     16   !!   Agrif_Update_Tra   : T-S agrif update 
     17   !!   Agrif_Update_Dyn   : dynamics agrif update 
     18   !!   Agrif_Update_ssh   : sea surface height update 
     19   !!   Agrif_Update_Tke   :  
     20   !!   Agrif_Update_vvl   :  
     21   !!   dom_vvl_update_UVF :  
     22   !!   updateTS           :  
     23   !!   updateu            : 
     24   !!   correct_u_bdy      : 
     25   !!   updatev            : 
     26   !!   correct_v_bdy      : 
     27   !!   updateu2d          : 
     28   !!   updatev2d          : 
     29   !!   updateSSH          : 
     30   !!   updateub2b         : 
     31   !!   reflux_sshu        : 
     32   !!   updatevb2b         : 
     33   !!   reflux_sshv        : 
     34   !!   update_scales      : 
     35   !!   updateEN           : 
     36   !!   updateAVT          : 
     37   !!   updateAVM          : 
     38   !!   updatee3t          : 
     39   !!---------------------------------------------------------------------- 
     40 
    1441#if defined key_agrif  
    1542   !!---------------------------------------------------------------------- 
    1643   !!   'key_agrif'                                              AGRIF zoom 
    1744   !!---------------------------------------------------------------------- 
    18    USE par_oce 
    19    USE oce 
    20    USE dom_oce 
     45   USE par_oce        ! ocean parameter 
     46   USE oce            ! ocean variables 
     47   USE dom_oce        ! ocean domain 
    2148   USE zdf_oce        ! vertical physics: ocean variables  
    22    USE agrif_oce 
     49   USE agrif_oce      !  
    2350   ! 
    2451   USE in_out_manager ! I/O manager 
    … …  
    6794      ! 
    6895   END SUBROUTINE Agrif_Update_Tra 
     96 
    6997 
    7098   SUBROUTINE Agrif_Update_Dyn( ) 
    … …  
    125153   END SUBROUTINE Agrif_Update_Dyn 
    126154 
     155 
    127156   SUBROUTINE Agrif_Update_ssh( ) 
    128       !!--------------------------------------------- 
    129       !!   *** ROUTINE Agrif_Update_ssh *** 
    130       !!--------------------------------------------- 
     157      !!---------------------------------------------------------------------- 
     158      !!                   *** ROUTINE Agrif_Update_ssh *** 
     159      !!---------------------------------------------------------------------- 
    131160      !  
    132161      IF (Agrif_Root()) RETURN 
    … …  
    163192 
    164193   SUBROUTINE Agrif_Update_Tke( ) 
    165       !!--------------------------------------------- 
    166       !!   *** ROUTINE Agrif_Update_Tke *** 
    167       !!--------------------------------------------- 
    168       !! 
     194      !!---------------------------------------------------------------------- 
     195      !!                   *** ROUTINE Agrif_Update_Tke *** 
     196      !!---------------------------------------------------------------------- 
    169197      !  
    170198      IF (Agrif_Root()) RETURN 
    171199      !        
    172200#  if defined TWO_WAY 
    173  
     201      ! 
    174202      Agrif_UseSpecialValueInUpdate = .TRUE. 
    175203      Agrif_SpecialValueFineGrid = 0. 
    176  
     204      ! 
    177205      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  ) 
    178206      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT ) 
    179207      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM ) 
    180  
     208      ! 
    181209      Agrif_UseSpecialValueInUpdate = .FALSE. 
    182  
     210      ! 
    183211#  endif 
    184        
     212      ! 
    185213   END SUBROUTINE Agrif_Update_Tke 
    186214 
    187215 
    188216   SUBROUTINE Agrif_Update_vvl( ) 
    189       !!--------------------------------------------- 
    190       !!   *** ROUTINE Agrif_Update_vvl *** 
    191       !!--------------------------------------------- 
    192       ! 
    193       IF (Agrif_Root()) RETURN 
     217      !!---------------------------------------------------------------------- 
     218      !!                   *** ROUTINE Agrif_Update_vvl *** 
     219      !!---------------------------------------------------------------------- 
     220      ! 
     221      IF ( Agrif_Root() )  RETURN 
    194222      ! 
    195223#if defined TWO_WAY   
    … …  
    214242   END SUBROUTINE Agrif_Update_vvl 
    215243 
     244 
    216245   SUBROUTINE dom_vvl_update_UVF 
    217       !!--------------------------------------------- 
    218       !!       *** ROUTINE dom_vvl_update_UVF *** 
    219       !!--------------------------------------------- 
    220       !! 
    221       INTEGER :: jk 
    222       REAL(wp):: zcoef 
    223       !!--------------------------------------------- 
    224  
    225       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 
    226                   & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
    227  
    228       ! Save "old" scale factor (prior update) for subsequent asselin correction 
    229       ! of prognostic variables 
     246      !!---------------------------------------------------------------------- 
     247      !!                   *** ROUTINE dom_vvl_update_UVF *** 
     248      !!---------------------------------------------------------------------- 
     249      INTEGER ::   jk      ! dummy loop index 
     250      REAL(wp)::   zcoef   ! local scalar 
     251      !!---------------------------------------------------------------------- 
     252      ! 
     253      IF (lwp.AND.lk_agrif_debug)   Write(*,*) 'Finalize e3 on grid Number', & 
     254         &                                      Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     255 
     256      ! Save "old" scale factor (prior update) for subsequent asselin correction of prognostic variables 
    230257      ! ----------------------- 
    231       ! 
    232258      e3u_a(:,:,:) = e3u_n(:,:,:) 
    233259      e3v_a(:,:,:) = e3v_n(:,:,:) 
    … …  
    239265      ! 1) NOW fields 
    240266      !-------------- 
    241        
    242          ! Vertical scale factor interpolations 
    243          ! ------------------------------------ 
    244       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' ) 
    245       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' ) 
    246       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' ) 
    247  
     267      !                       ! Vertical scale factor interpolations 
     268      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n (:,:,:) , 'U' ) 
     269      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n (:,:,:) , 'V' ) 
     270      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n (:,:,:) , 'F' ) 
    248271      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    249272      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    250  
    251          ! Update total depths: 
    252          ! -------------------- 
     273      ! 
     274      !                       ! Update total depths 
    253275      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points 
    254276      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points 
    … …  
    264286      ! 2) BEFORE fields: 
    265287      !------------------ 
    266       IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
    267          ! 
    268          ! Vertical scale factor interpolations 
    269          ! ------------------------------------ 
    270          CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  ) 
    271          CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  ) 
    272  
     288      IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     289!!gm      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     290         !                    ! Vertical scale factor interpolations 
     291         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b (:,:,:), 'U'  ) 
     292         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b (:,:,:), 'V'  ) 
    273293         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    274294         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    275  
    276          ! Update total depths: 
    277          ! -------------------- 
     295         ! 
     296         !                    ! Update total depths: 
    278297         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points 
    279298         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points 
    … …  
    289308   END SUBROUTINE dom_vvl_update_UVF 
    290309 
    291 #if defined key_vertical 
     310# if defined key_vertical 
    292311 
    293312   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    294313      !!---------------------------------------------------------------------- 
    295       !!           *** ROUTINE updateT *** 
    296       !!--------------------------------------------- 
     314      !!                   *** ROUTINE updateT *** 
     315      !!---------------------------------------------------------------------- 
    297316      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    298317      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    306325      REAL(wp) :: zrho_xy, h_diff 
    307326      REAL(wp) :: tabin(k1:k2,n1:n2) 
    308       !!--------------------------------------------- 
     327      !!---------------------------------------------------------------------- 
    309328      ! 
    310329      IF (before) THEN 
    311330         AGRIF_SpecialValue = -999._wp 
    312331         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
    313          DO jn = n1,n2-1 
    314             DO jk=k1,k2 
    315                DO jj=j1,j2 
    316                   DO ji=i1,i2 
     332         DO jn = n1, n2-1 
     333            DO jk = k1, k2 
     334               DO jj = j1, j2 
     335                  DO ji = i1, i2 
    317336                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
    318337                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
    … …  
    321340            END DO 
    322341         END DO 
    323          DO jk=k1,k2 
    324             DO jj=j1,j2 
    325                DO ji=i1,i2 
     342         DO jk = k1, k2 
     343            DO jj = j1, j2 
     344               DO ji = i1, i2 
    326345                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
    327346                                           + (tmask(ji,jj,jk)-1)*999._wp 
    … …  
    332351         tabres_child(:,:,:,:) = 0. 
    333352         AGRIF_SpecialValue = 0._wp 
    334          DO jj=j1,j2 
    335             DO ji=i1,i2 
     353         DO jj = j1 , j2 
     354            DO ji = i1, i2 
    336355               N_in = 0 
    337                DO jk=k1,k2 !k2 = jpk of child grid 
    338                   IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
     356               DO jk = k1, k2 !k2 = jpk of child grid 
     357                  IF ( tabres(ji,jj,jk,n2) == 0  )  EXIT 
    339358                  N_in = N_in + 1 
    340359                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    341                   h_in(N_in) = tabres(ji,jj,jk,n2) 
    342                ENDDO 
     360                  h_in (N_in) = tabres(ji,jj,jk,n2) 
     361               END DO 
    343362               N_out = 0 
    344                DO jk=1,jpk ! jpk of parent grid 
    345                   IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
     363               DO jk = 1, jpk ! jpk of parent grid 
     364                  IF (tmask(ji,jj,jk) < -900)   EXIT ! TODO: Will not work with ISF 
    346365                  N_out = N_out + 1 
    347366                  h_out(N_out) = e3t_n(ji,jj,jk)  
    348                ENDDO 
     367               END DO 
    349368               IF (N_in > 0) THEN !Remove this? 
    350369                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    … …  
    355374                     STOP 
    356375                  ENDIF 
    357                   DO jn=n1,n2-1 
     376                  DO jn = n1, n2-1 
    358377                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    359                   ENDDO 
     378                  END DO 
    360379               ENDIF 
    361             ENDDO 
    362          ENDDO 
    363  
    364          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     380            END DO 
     381         END DO 
     382 
     383         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     384!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    365385            ! Add asselin part 
    366             DO jn = n1,n2-1 
    367                DO jk=1,jpk 
    368                   DO jj=j1,j2 
    369                      DO ji=i1,i2 
    370                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
     386            DO jn = n1, n2-1 
     387               DO jk = 1, jpk 
     388                  DO jj = j1, j2 
     389                     DO ji = i1, i2 
     390                        IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN 
    371391                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    372392                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    373393                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    374394                        ENDIF 
    375                      ENDDO 
    376                   ENDDO 
    377                ENDDO 
    378             ENDDO 
    379          ENDIF 
    380          DO jn = n1,n2-1 
    381             DO jk=1,jpk 
    382                DO jj=j1,j2 
    383                   DO ji=i1,i2 
    384                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     395                     END DO 
     396                  END DO 
     397               END DO 
     398            END DO 
     399         ENDIF 
     400         DO jn = n1, n2-1 
     401            DO jk = 1, jpk 
     402               DO jj = j1, j2 
     403                  DO ji = i1, i2 
     404                     IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN  
    385405                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    386406                     END IF 
    … …  
    396416 
    397417   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    398       !!--------------------------------------------- 
    399       !!           *** ROUTINE updateT *** 
    400       !!--------------------------------------------- 
     418      !!---------------------------------------------------------------------- 
     419      !!                   *** ROUTINE ROUTINE updateT *** 
     420      !!---------------------------------------------------------------------- 
    401421      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    402422      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    403423      LOGICAL, INTENT(in) :: before 
    404       !! 
     424      ! 
    405425      INTEGER :: ji,jj,jk,jn 
    406426      REAL(wp) :: ztb, ztnu, ztno 
    407       !!--------------------------------------------- 
     427      !!---------------------------------------------------------------------- 
    408428      ! 
    409429      IF (before) THEN 
    … …  
    425445            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    426446                                         & * tmask(i1:i2,j1:j2,k1:k2) 
    427          ENDDO 
     447         END DO 
    428448!< jc tmp 
    429          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     449         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     450!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    430451            ! Add asselin part 
    431452            DO jn = 1,jpts 
    … …  
    457478         END DO 
    458479         ! 
    459          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     480         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     481!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    460482            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts) 
    461483         ENDIF 
    … …  
    470492 
    471493   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    472       !!--------------------------------------------- 
    473       !!           *** ROUTINE updateu *** 
    474       !!--------------------------------------------- 
     494      !!---------------------------------------------------------------------- 
     495      !!                   *** ROUTINE updateu *** 
     496      !!---------------------------------------------------------------------- 
    475497      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    476498      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    487509      REAL(wp) :: tabin(k1:k2) 
    488510! VERTICAL REFINEMENT END 
    489       !!--------------------------------------------- 
     511      !!---------------------------------------------------------------------- 
    490512      !  
    491513      IF( before ) THEN 
    … …  
    515537                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    516538                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
    517                ENDDO 
     539               END DO 
    518540               N_out = 0 
    519541               DO jk=1,jpk 
    … …  
    521543                  N_out = N_out + 1 
    522544                  h_out(N_out) = e3u_n(ji,jj,jk) 
    523                ENDDO 
     545               END DO 
    524546               IF (N_in * N_out > 0) THEN 
    525547                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    … …  
    538560                           EXIT 
    539561                        ENDIF 
    540                      ENDDO 
     562                     END DO 
    541563                  ENDIF 
    542564                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    543565                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
    544566               ENDIF 
    545             ENDDO 
    546          ENDDO 
     567            END DO 
     568         END DO 
    547569 
    548570         DO jk=1,jpk 
    549571            DO jj=j1,j2 
    550572               DO ji=i1,i2 
    551                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     573                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler) ) THEN    ! Add asselin part 
     574!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
    552575                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
    553576                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
    … …  
    565588 
    566589   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    567       !!--------------------------------------------- 
    568       !!           *** ROUTINE updateu *** 
    569       !!--------------------------------------------- 
     590      !!---------------------------------------------------------------------- 
     591      !!                   *** ROUTINE updateu *** 
     592      !!---------------------------------------------------------------------- 
    570593      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    571594      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    574597      INTEGER  :: ji, jj, jk 
    575598      REAL(wp) :: zrhoy, zub, zunu, zuno 
    576       !!--------------------------------------------- 
     599      !!---------------------------------------------------------------------- 
    577600      !  
    578601      IF( before ) THEN 
    … …  
    587610                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    588611                  ! 
    589                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     612                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     613!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
    590614                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used 
    591615                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 
    … …  
    600624         END DO 
    601625         ! 
    602          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     626         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     627!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    603628            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2) 
    604629         ENDIF 
    … …  
    611636 
    612637   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    613       !!--------------------------------------------- 
    614       !!           *** ROUTINE correct_u_bdy *** 
    615       !!--------------------------------------------- 
     638      !!---------------------------------------------------------------------- 
     639      !!                   *** ROUTINE correct_u_bdy *** 
     640      !!---------------------------------------------------------------------- 
    616641      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    617642      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    618643      LOGICAL                                     , INTENT(in   ) :: before 
    619       INTEGER                                     , INTENT(in)    :: nb, ndir 
     644      INTEGER                                     , INTENT(in   ) :: nb, ndir 
    620645      !! 
    621646      LOGICAL :: western_side, eastern_side  
    622       ! 
    623       INTEGER  :: jj, jk 
    624       REAL(wp) :: zcor 
    625       !!--------------------------------------------- 
     647      INTEGER ::   jj, jk 
     648      REAL(wp)::   zcor 
     649      !!---------------------------------------------------------------------- 
    626650      !  
    627651      IF( .NOT.before ) THEN 
    … …  
    657681 
    658682   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    659       !!--------------------------------------------- 
    660       !!           *** ROUTINE updatev *** 
    661       !!--------------------------------------------- 
     683      !!---------------------------------------------------------------------- 
     684      !!                   *** ROUTINE updatev *** 
     685      !!---------------------------------------------------------------------- 
    662686      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    663687      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    674698      REAL(wp) :: tabin(k1:k2) 
    675699! VERTICAL REFINEMENT END 
    676       !!---------------------------------------------       
     700      !!----------------------------------------------------------------------       
    677701      ! 
    678702      IF( before ) THEN 
    … …  
    700724                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    701725                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
    702                ENDDO 
     726               END DO 
    703727               N_out = 0 
    704728               DO jk=1,jpk 
    … …  
    706730                  N_out = N_out + 1 
    707731                  h_out(N_out) = e3v_n(ji,jj,jk) 
    708                ENDDO 
     732               END DO 
    709733               IF (N_in * N_out > 0) THEN 
    710734                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    … …  
    723747                           EXIT 
    724748                        ENDIF 
    725                      ENDDO 
     749                     END DO 
    726750                  ENDIF 
    727751                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    728752                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    729753               ENDIF 
    730             ENDDO 
    731          ENDDO 
     754            END DO 
     755         END DO 
    732756 
    733757         DO jk=1,jpk 
    … …  
    735759               DO ji=i1,i2 
    736760                  ! 
    737                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 
     761                  IF( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN  ! Add asselin part 
     762!!gm                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN  ! Add asselin part 
    738763                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    739764                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    … …  
    751776 
    752777   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
    753       !!--------------------------------------------- 
    754       !!           *** ROUTINE updatev *** 
    755       !!--------------------------------------------- 
     778      !!---------------------------------------------------------------------- 
     779      !!                   *** ROUTINE updatev *** 
     780      !!---------------------------------------------------------------------- 
    756781      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    757782      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    760785      INTEGER  :: ji, jj, jk 
    761786      REAL(wp) :: zrhox, zvb, zvnu, zvno 
    762       !!---------------------------------------------       
     787      !!----------------------------------------------------------------------       
    763788      ! 
    764789      IF (before) THEN 
    … …  
    777802                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
    778803                  ! 
    779                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     804                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     805!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
    780806                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 
    781807                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 
    … …  
    790816         END DO 
    791817         ! 
    792          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     818         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     819!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    793820            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2) 
    794821         ENDIF 
    … …  
    801828 
    802829   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    803       !!--------------------------------------------- 
    804       !!           *** ROUTINE correct_u_bdy *** 
    805       !!--------------------------------------------- 
     830      !!---------------------------------------------------------------------- 
     831      !!                   *** ROUTINE correct_v_bdy *** 
     832      !!---------------------------------------------------------------------- 
    806833      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    807834      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    … …  
    813840      INTEGER  :: ji, jk 
    814841      REAL(wp) :: zcor 
    815       !!--------------------------------------------- 
     842      !!---------------------------------------------------------------------- 
    816843      !  
    817844      IF( .NOT.before ) THEN 
    … …  
    847874   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before ) 
    848875      !!---------------------------------------------------------------------- 
    849       !!                      *** ROUTINE updateu2d *** 
     876      !!                   *** ROUTINE updateu2d *** 
    850877      !!---------------------------------------------------------------------- 
    851878      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    852879      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    853880      LOGICAL                         , INTENT(in   ) ::   before 
    854       !!  
     881      ! 
    855882      INTEGER  :: ji, jj, jk 
    856883      REAL(wp) :: zrhoy 
    857884      REAL(wp) :: zcorr 
    858       !!--------------------------------------------- 
     885      !!---------------------------------------------------------------------- 
    859886      ! 
    860887      IF( before ) THEN 
    … …  
    883910               ! Update barotropic velocities: 
    884911               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    885                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     912                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     913!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    886914                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 
    887915                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
    … …  
    904932         END DO 
    905933         ! 
    906          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     934         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     935!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    907936            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2) 
    908937         ENDIF 
    … …  
    948977               ! 
    949978               ! Update barotropic velocities: 
    950                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    951                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     979               IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN 
     980                  IF ( .NOT.( lk_agrif_fstep. AND. l_1st_euler ) ) THEN    ! Add asselin part 
     981!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    952982                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 
    953983                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
    … …  
    9701000         END DO 
    9711001         ! 
    972          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1002         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     1003!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    9731004            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2) 
    9741005         ENDIF 
    … …  
    9861017      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    9871018      LOGICAL                         , INTENT(in   ) ::   before 
    988       !! 
     1019      ! 
    9891020      INTEGER :: ji, jj 
    9901021      !!---------------------------------------------------------------------- 
    … …  
    9971028         END DO 
    9981029      ELSE 
    999          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     1030         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     1031!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    10001032            DO jj=j1,j2 
    10011033               DO ji=i1,i2 
    … …  
    10121044         END DO 
    10131045         ! 
    1014          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1046         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     1047!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    10151048            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2) 
    10161049         ENDIF 
    10171050         ! 
    1018  
    10191051      ENDIF 
    10201052      ! 
    … …  
    10621094   END SUBROUTINE updateub2b 
    10631095 
     1096 
    10641097   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir ) 
    1065       !!--------------------------------------------- 
    1066       !!          *** ROUTINE reflux_sshu *** 
    1067       !!--------------------------------------------- 
    1068       INTEGER, INTENT(in) :: i1, i2, j1, j2 
    1069       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
    1070       LOGICAL, INTENT(in) :: before 
    1071       INTEGER, INTENT(in) :: nb, ndir 
    1072       !! 
    1073       LOGICAL :: western_side, eastern_side  
    1074       INTEGER :: ji, jj 
    1075       REAL(wp) :: zrhoy, za1, zcor 
    1076       !!--------------------------------------------- 
     1098      !!---------------------------------------------------------------------- 
     1099      !!                   *** ROUTINE reflux_sshu *** 
     1100      !!---------------------------------------------------------------------- 
     1101      INTEGER                         , INTENT(in   ) ::  i1, i2, j1, j2 
     1102      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
     1103      LOGICAL                         , INTENT(in   ) ::  before 
     1104      INTEGER                         , INTENT(in   ) ::  nb, ndir 
     1105      ! 
     1106      LOGICAL ::   western_side, eastern_side  
     1107      INTEGER ::   ji, jj 
     1108      REAL(wp)::  zrhoy, za1, zcor 
     1109      !!---------------------------------------------------------------------- 
    10771110      ! 
    10781111      IF (before) THEN 
    … …  
    10911124         eastern_side  = (nb == 1).AND.(ndir == 2) 
    10921125         ! 
    1093          IF (western_side) THEN 
     1126         IF ( western_side ) THEN 
    10941127            DO jj=j1,j2 
    10951128               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))  
    10961129               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor 
    1097                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
    1098             END DO 
    1099          ENDIF 
    1100          IF (eastern_side) THEN 
     1130               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
     1131!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
     1132            END DO 
     1133         ENDIF 
     1134         IF ( eastern_side ) THEN 
    11011135            DO jj=j1,j2 
    11021136               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 
    11031137               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor 
    1104                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
     1138               IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
     1139!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
    11051140            END DO 
    11061141         ENDIF 
    … …  
    11101145   END SUBROUTINE reflux_sshu 
    11111146 
     1147 
    11121148   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before ) 
    11131149      !!---------------------------------------------------------------------- 
    1114       !!                      *** ROUTINE updatevb2b *** 
     1150      !!                    *** ROUTINE updatevb2b *** 
    11151151      !!---------------------------------------------------------------------- 
    11161152      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    11171153      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    11181154      LOGICAL                         , INTENT(in   ) ::   before 
    1119       !! 
     1155      ! 
    11201156      INTEGER :: ji, jj 
    11211157      REAL(wp) :: zrhox, za1, zcor 
    1122       !!--------------------------------------------- 
     1158      !!--------------------------------------------------------------------- 
    11231159      ! 
    11241160      IF( before ) THEN 
    … …  
    11501186   END SUBROUTINE updatevb2b 
    11511187 
     1188 
    11521189   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir ) 
    1153       !!--------------------------------------------- 
    1154       !!          *** ROUTINE reflux_sshv *** 
    1155       !!--------------------------------------------- 
    1156       INTEGER, INTENT(in) :: i1, i2, j1, j2 
    1157       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
    1158       LOGICAL, INTENT(in) :: before 
    1159       INTEGER, INTENT(in) :: nb, ndir 
     1190      !!---------------------------------------------------------------------- 
     1191      !!                   *** ROUTINE reflux_sshv *** 
     1192      !!---------------------------------------------------------------------- 
     1193      INTEGER                         , INTENT(in   ) ::  i1, i2, j1, j2 
     1194      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
     1195      LOGICAL                         , INTENT(in   ) ::  before 
     1196      INTEGER                         , INTENT(in   ) ::  nb, ndir 
    11601197      !! 
    11611198      LOGICAL :: southern_side, northern_side  
    11621199      INTEGER :: ji, jj 
    11631200      REAL(wp) :: zrhox, za1, zcor 
    1164       !!--------------------------------------------- 
     1201      !!---------------------------------------------------------------------- 
    11651202      ! 
    11661203      IF (before) THEN 
    … …  
    11831220               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1)) 
    11841221               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor 
    1185                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
     1222               IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) )   sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
     1223!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
    11861224            END DO 
    11871225         ENDIF 
    … …  
    11901228               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2)) 
    11911229               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor 
    1192                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
     1230               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
     1231!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
    11931232            END DO 
    11941233         ENDIF 
    … …  
    11981237   END SUBROUTINE reflux_sshv 
    11991238 
     1239 
    12001240   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before ) 
     1241      !!---------------------------------------------------------------------- 
     1242      !!                      *** ROUTINE updateT *** 
    12011243      ! 
    12021244      ! ====>>>>>>>>>>    currently not used 
    12031245      ! 
    1204       !!---------------------------------------------------------------------- 
    1205       !!                      *** ROUTINE updateT *** 
    12061246      !!---------------------------------------------------------------------- 
    12071247      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    … …  
    12841324 
    12851325   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before ) 
    1286       !!--------------------------------------------- 
    1287       !!           *** ROUTINE updateavm *** 
     1326      !!---------------------------------------------------------------------- 
     1327      !!                      *** ROUTINE updateavm *** 
    12881328      !!---------------------------------------------------------------------- 
    12891329      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    … …  
    12981338   END SUBROUTINE updateAVM 
    12991339 
     1340 
    13001341   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) 
    1301       !!--------------------------------------------- 
    1302       !!           *** ROUTINE updatee3t *** 
    1303       !!--------------------------------------------- 
     1342      !!---------------------------------------------------------------------- 
     1343      !!                   *** ROUTINE updatee3t *** 
     1344      !!---------------------------------------------------------------------- 
    13041345      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum 
    13051346      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 
    … …  
    13131354      IF (.NOT.before) THEN 
    13141355         ! 
    1315          ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))  
     1356         ALLOCATE( ptab(i1:i2,j1:j2,1:jpk) )  
    13161357         ! 
    13171358         ! Update e3t from ssh (z* case only) 
    … …  
    13351376!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1) 
    13361377 
    1337          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     1378         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler==0 ) ) THEN 
     1379!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
    13381380            DO jk = 1, jpkm1 
    13391381               DO jj=j1,j2 
    … …  
    13981440         END DO 
    13991441         ! 
    1400          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1442         IF ( l_1st_euler .AND. Agrif_Nb_Step()==0 ) THEN 
     1443!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    14011444            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk) 
    14021445            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/NST/agrif_top_update.F90

    r9598 r9863  
    109109                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    110110                  h_in(N_in) = tabres(ji,jj,jk,n2) 
    111                ENDDO 
     111               END DO 
    112112               N_out = 0 
    113113               DO jk=1,jpk ! jpk of parent grid 
    … …  
    115115                  N_out = N_out + 1 
    116116                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
    117                ENDDO 
     117               END DO 
    118118               IF (N_in > 0) THEN !Remove this? 
    119119                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    … …  
    126126                  DO jn=1,jptra 
    127127                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    128                   ENDDO 
     128                  END DO 
    129129               ENDIF 
    130             ENDDO 
    131          ENDDO 
    132  
    133          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     130            END DO 
     131         END DO 
     132 
     133         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     134!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    134135            ! Add asselin part 
    135136            DO jn = 1,jptra 
    … …  
    142143                                 &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    143144                        ENDIF 
    144                      ENDDO 
    145                   ENDDO 
    146                ENDDO 
    147             ENDDO 
     145                     END DO 
     146                  END DO 
     147               END DO 
     148            END DO 
    148149         ENDIF 
    149150         DO jn = 1,jptra 
    … …  
    195196            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    196197                                         & * tmask(i1:i2,j1:j2,k1:k2) 
    197          ENDDO 
     198         END DO 
    198199!< jc tmp 
    199          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     200         IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     201!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    200202            ! Add asselin part 
    201203            DO jn = n1,n2 
    … …  
    210212                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
    211213                        ENDIF 
    212                      ENDDO 
    213                   ENDDO 
    214                ENDDO 
    215             ENDDO 
     214                     END DO 
     215                  END DO 
     216               END DO 
     217            END DO 
    216218         ENDIF 
    217219         DO jn = n1,n2 
    … …  
    227229         END DO 
    228230         ! 
    229          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     231         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     232!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    230233            trb(i1:i2,j1:j2,k1:k2,n1:n2)  = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
    231234         ENDIF 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/ASM/asminc.F90

    r9656 r9863  
    491491      ENDIF 
    492492      ! 
    493       IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler 
     493      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', ln_1st_euler 
    494494      ! 
    495495      IF( lk_asminc ) THEN                            !==  data assimilation  ==! 
    … …  
    579579         IF ( kt == nitdin_r ) THEN 
    580580            ! 
    581             neuler = 0  ! Force Euler forward step 
     581            l_1st_euler = .TRUE.  ! Force Euler forward step 
    582582            ! 
    583583            ! Initialize the now fields with the background + increment 
    … …  
    677677         IF ( kt == nitdin_r ) THEN 
    678678            ! 
    679             neuler = 0                    ! Force Euler forward step 
     679            l_1st_euler = .TRUE.          ! Force Euler forward step 
    680680            ! 
    681681            ! Initialize the now fields with the background + increment 
    … …  
    752752         IF ( kt == nitdin_r ) THEN 
    753753            ! 
    754             neuler = 0                                   ! Force Euler forward step 
     754            l_1st_euler = .TRUE.                         ! Force Euler forward step 
    755755            ! 
    756756            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment 
    … …  
    758758            sshb(:,:) = sshn(:,:)                        ! Update before fields 
    759759            e3t_b(:,:,:) = e3t_n(:,:,:) 
    760 !!gm why not e3u_b, e3v_b, gdept_b ???? 
     760             
     761!!gm BUG :   missing the update of all other scale factors (e3u e3v e3w  etc... _n and _b)  
     762!!           see dom_vvl_init  
    761763            ! 
    762764            DEALLOCATE( ssh_bkg    ) 
    … …  
    894896         IF ( kt == nitdin_r ) THEN 
    895897            ! 
    896             neuler = 0                    ! Force Euler forward step 
     898            l_1st_euler = .TRUE.             ! Force Euler forward step 
    897899            ! 
    898900            ! Sea-ice : SI3 case 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diacfl.F90

    r9598 r9863  
    5555      ! 
    5656      INTEGER :: ji, jj, jk   ! dummy loop indices 
    57       REAL(wp)::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
     57      REAL(wp)::   zCu_max, zCv_max, zCw_max       ! local scalars 
    5858      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc   ! workspace 
    5959!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
    … …  
    6262      IF( ln_timing )   CALL timing_start('dia_cfl') 
    6363      ! 
    64       !                       ! setup timestep multiplier to account for initial Eulerian timestep 
    65       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2dt = rdt 
    66       ELSE                                        ;    z2dt = rdt * 2._wp 
    67       ENDIF 
    68       ! 
    6964      !                 
    7065      DO jk = 1, jpk       ! calculate Courant numbers 
    7166         DO jj = 1, jpj 
    7267            DO ji = 1, fs_jpim1   ! vector opt. 
    73                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    74                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
    75                zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction 
     68               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * r2dt / e1u  (ji,jj)      ! for i-direction 
     69               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * r2dt / e2v  (ji,jj)      ! for j-direction 
     70               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * r2dt / e3w_n(ji,jj,jk)   ! for k-direction 
    7671            END DO 
    7772         END DO          
    … …  
    120115         WRITE(numcfl,*) '******************************************' 
    121116         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3) 
    122          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max 
     117         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCu_max 
    123118         WRITE(numcfl,*) '******************************************' 
    124119         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3) 
    125          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max 
     120         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCv_max 
    126121         WRITE(numcfl,*) '******************************************' 
    127122         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3) 
    128          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max 
     123         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCw_max 
    129124         CLOSE( numcfl )  
    130125         ! 
    … …  
    133128         WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run ' 
    134129         WRITE(numout,*) '~~~~~~~' 
    135          WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max 
    136          WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max 
    137          WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max 
     130         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', r2dt/rCu_max 
     131         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', r2dt/rCv_max 
     132         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', r2dt/rCw_max 
    138133      ENDIF 
    139134      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/dom_oce.F90

    r9667 r9863  
    3535   REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer 
    3636   REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter 
    37    INTEGER , PUBLIC ::   nn_euler       !: =0 start with forward time step or not (=1) 
     37   LOGICAL , PUBLIC ::   ln_1st_euler   !: =0 start with forward time step or not (=1) 
    3838   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    3939   LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers 
    … …  
    5757   !                                   !! old non-DOCTOR names still used in the model 
    5858   REAL(wp), PUBLIC ::   atfp           !: asselin time filter parameter 
    59    REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer 
     59   REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer (=rn_rdt) 
    6060 
    6161   !                                   !!! associated variables 
    62    INTEGER , PUBLIC ::   neuler         !: restart euler forward option (0=Euler) 
    63    REAL(wp), PUBLIC ::   r2dt           !: = 2*rdt except at nit000 (=rdt) if neuler=0 
     62   LOGICAL , PUBLIC ::   l_1st_euler    !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T) 
     63   REAL(wp), PUBLIC ::   r2dt, r1_2dt   !: = 2*rdt and 1/(2*rdt) except if l_1st_euler=T) 
    6464 
    6565   !!---------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domain.F90

    r9598 r9863  
    288288      INTEGER  ::   ios   ! Local integer 
    289289      ! 
    290       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
    291          &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     & 
    292          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
    293          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
     290      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                   & 
     291         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl   ,     & 
     292         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate   ,     & 
     293         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler,     & 
    294294         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
    295295      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
    … …  
    323323         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir ) 
    324324         WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart 
    325          WRITE(numout,*) '      start with forward time step    nn_euler        = ', nn_euler 
     325         WRITE(numout,*) '      start with forward time step    ln_1st_euler    = ', ln_1st_euler 
    326326         WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl 
    327327         WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000 
    … …  
    361361      nstocklist = nn_stocklist 
    362362      nwrite = nn_write 
    363       neuler = nn_euler 
    364       IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN 
     363      IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward)  
     364         l_1st_euler = ln_1st_euler 
     365      ELSE                       ! start from rest : always an Euler scheme for the 1st time-step 
    365366         IF(lwp) WRITE(numout,*)   
    366367         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
    367          IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : nn_euler is forced to 0 '    
    368          neuler = 0 
     368         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '    
     369         l_1st_euler = .TRUE. 
    369370      ENDIF 
    370371      !                             ! control of output frequency 
    … …  
    374375         nstock = nitend 
    375376      ENDIF 
    376       IF ( nwrite == 0 ) THEN 
     377      IF( nwrite == 0 ) THEN 
    377378         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend 
    378379         CALL ctl_warn( ctmp1 ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90

    r9598 r9863  
    5656   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
    5757   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence 
    58    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors 
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors 
     58   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors 
     59   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a   ! baroclinic scale factors 
    6060   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors 
    6161   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    … …  
    7676      IF( ln_vvl_zstar )   dom_vvl_alloc = 0 
    7777      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    78          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
    79             &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
     78         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   & 
     79            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
    8080            &      STAT = dom_vvl_alloc        ) 
    8181         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    … …  
    103103      !!               - interpolate scale factors 
    104104      !! 
    105       !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     105      !! ** Action  : - e3t_(n/b) and te3t_(n/b) 
    106106      !!              - Regrid: e3(u/v)_n 
    107107      !!                        e3(u/v)_b        
    … …  
    129129      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    130130      ! 
    131       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     131      !                    ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf 
    132132      CALL dom_vvl_rst( nit000, 'READ' ) 
    133133      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    … …  
    280280      !! 
    281281      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
    282       !!               - tilde_e3t_a: after increment of vertical scale factor  
     282      !!               - te3t_a: after increment of vertical scale factor  
    283283      !!                              in z_tilde case 
    284284      !!               - e3(t/u/v)_a 
    … …  
    353353         ! II - after z_tilde increments of vertical scale factors 
    354354         ! ======================================================= 
    355          tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms 
     355         te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms 
    356356 
    357357         ! 1 - High frequency divergence term 
    … …  
    359359         IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    360360            DO jk = 1, jpkm1 
    361                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     361               te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    362362            END DO 
    363363         ELSE                         ! layer case 
    364364            DO jk = 1, jpkm1 
    365                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     365               te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    366366            END DO 
    367367         ENDIF 
    … …  
    371371         IF( ln_vvl_ztilde ) THEN 
    372372            DO jk = 1, jpk 
    373                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
     373               te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 
    374374            END DO 
    375375         ENDIF 
    … …  
    383383               DO ji = 1, fs_jpim1   ! vector opt. 
    384384                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    385                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     385                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) ) 
    386386                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    387                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     387                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji  ,jj+1,jk) ) 
    388388                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    389389                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    … …  
    400400            DO jj = 2, jpjm1 
    401401               DO ji = fs_2, fs_jpim1   ! vector opt. 
    402                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     402                  te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    403403                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    404404                     &                                            ) * r1_e1e2t(ji,jj) 
    … …  
    414414         ! Leapfrog time stepping 
    415415         ! ~~~~~~~~~~~~~~~~~~~~~~ 
    416          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    417             z2dt =  rdt 
    418          ELSE 
    419             z2dt = 2.0_wp * rdt 
    420          ENDIF 
    421          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    422          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     416         CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 
     417         te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 
    423418 
    424419         ! Maximum deformation control 
    … …  
    426421         ze3t(:,:,jpk) = 0._wp 
    427422         DO jk = 1, jpkm1 
    428             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
     423            ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
    429424         END DO 
    430425         z_tmax = MAXVAL( ze3t(:,:,:) ) 
    … …  
    446441            ENDIF 
    447442            IF (lwp) THEN 
    448                WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
     443               WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
    449444               WRITE(numout, *) 'at i, j, k=', ijk_max 
    450                WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
     445               WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
    451446               WRITE(numout, *) 'at i, j, k=', ijk_min             
    452                CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
     447               CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
    453448            ENDIF 
    454449         ENDIF 
    455450         ! - ML - end test 
    456451         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    457          tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
    458          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
     452         te3t_a(:,:,:) = MIN( te3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
     453         te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    459454 
    460455         ! 
    … …  
    462457         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    463458         DO jk = 1, jpkm1 
    464             dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 
     459            dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 
    465460         END DO 
    466461         ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    … …  
    470465         !        i.e. locally and not spread over the water column. 
    471466         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 
    472          zht(:,:) = 0. 
     467         zht(:,:) = 0._wp 
    473468         DO jk = 1, jpkm1 
    474             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     469            zht(:,:)  = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 
    475470         END DO 
    476471         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    477472         DO jk = 1, jpkm1 
    478             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     473            dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    479474         END DO 
    480475 
    … …  
    484479      !                                           ! ---baroclinic part--------- ! 
    485480         DO jk = 1, jpkm1 
    486             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     481            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 
    487482         END DO 
    488483      ENDIF 
    … …  
    494489            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 
    495490            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
    496             IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 
     491            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 
    497492         END IF 
    498493         ! 
    … …  
    573568      !!               - recompute depths and water height fields 
    574569      !! 
    575       !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 
     570      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step 
    576571      !!               - Recompute: 
    577572      !!                    e3(u/v)_b        
    … …  
    587582      INTEGER, INTENT( in ) ::   kt   ! time step 
    588583      ! 
    589       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    590       REAL(wp) ::   zcoef        ! local scalar 
     584      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     585      REAL(wp) ::   zcoef, ze3f   ! local scalar 
    591586      !!---------------------------------------------------------------------- 
    592587      ! 
    … …  
    605600      ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 
    606601      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    607          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    608             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 
     602         IF( l_1st_euler ) THEN 
     603            te3t_n(:,:,:) = te3t_a(:,:,:) 
    609604         ELSE 
    610             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) &  
    611             &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 
     605            DO jk = 1, jpk 
     606               DO jj = 1, jpj 
     607                  DO ji = 1, jpi 
     608                     ze3f = te3t_n(ji,jj,jk)   & 
     609                        & + atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 
     610                     te3t_b(ji,jj,jk) = ze3f 
     611                     te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 
     612                  END DO 
     613               END DO 
     614            END DO 
    612615         ENDIF 
    613          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    614616      ENDIF 
    615617      gdept_b(:,:,:) = gdept_n(:,:,:) 
    … …  
    806808            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
    807809            ! 
    808             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
    809             id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. ) 
     810            id1 = iom_varid( numror, 'e3t_b'      , ldstop = .FALSE. ) 
     811            id2 = iom_varid( numror, 'e3t_n'      , ldstop = .FALSE. ) 
    810812            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    811813            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    812             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     814            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. ) 
    813815            !                             ! --------- ! 
    814816            !                             ! all cases ! 
    … …  
    823825                  e3t_b(:,:,:) = e3t_0(:,:,:) 
    824826               END WHERE 
    825                IF( neuler == 0 ) THEN 
     827               IF( l_1st_euler ) THEN 
    826828                  e3t_b(:,:,:) = e3t_n(:,:,:) 
    827829               ENDIF 
    … …  
    829831               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
    830832               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    831                IF(lwp) write(numout,*) 'neuler is forced to 0' 
     833               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    832834               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    833835               e3t_n(:,:,:) = e3t_b(:,:,:) 
    834                neuler = 0 
     836               l_1st_euler = .TRUE. 
    835837            ELSE IF( id2 > 0 ) THEN 
    836838               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
    837839               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    838                IF(lwp) write(numout,*) 'neuler is forced to 0' 
     840               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    839841               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
    840842               e3t_b(:,:,:) = e3t_n(:,:,:) 
    841                neuler = 0 
     843               l_1st_euler = .TRUE. 
    842844            ELSE 
    843845               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
    844846               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    845                IF(lwp) write(numout,*) 'neuler is forced to 0' 
     847               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    846848               DO jk = 1, jpk 
    847849                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    … …  
    850852               END DO 
    851853               e3t_b(:,:,:) = e3t_n(:,:,:) 
    852                neuler = 0 
     854               l_1st_euler = .TRUE. 
    853855            ENDIF 
    854856            !                             ! ----------- ! 
    … …  
    862864               !                          ! ----------------------- ! 
    863865               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
    864                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios ) 
    865                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios ) 
     866                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios ) 
     867                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios ) 
    866868               ELSE                            ! one at least array is missing 
    867                   tilde_e3t_b(:,:,:) = 0.0_wp 
    868                   tilde_e3t_n(:,:,:) = 0.0_wp 
     869                  te3t_b(:,:,:) = 0.0_wp 
     870                  te3t_n(:,:,:) = 0.0_wp 
    869871               ENDIF 
    870872               !                          ! ------------ ! 
    … …  
    942944 
    943945            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    944                tilde_e3t_b(:,:,:) = 0._wp 
    945                tilde_e3t_n(:,:,:) = 0._wp 
     946               te3t_b(:,:,:) = 0._wp 
     947               te3t_n(:,:,:) = 0._wp 
    946948               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 
    947949            END IF 
    … …  
    960962         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
    961963            !                                        ! ----------------------- ! 
    962             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios) 
    963             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios) 
     964            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 
     965            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 
    964966         END IF 
    965967         !                                           ! -------------!     

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/iscplrst.F90

    r9598 r9863  
    8989      END IF 
    9090      ! 
    91       neuler = 0              ! next step is an euler time step 
     91      l_1st_euler = .TRUE.    ! next step is an euler time step 
    9292      ! 
    9393      !                       ! set _b and _n variables equal 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/istate.F90

    r9598 r9863  
    9292         !                                    ! --------------- 
    9393         numror = 0                           ! define numror = 0 -> no restart file to read 
    94          neuler = 0                           ! Set time-step indicator at nit000 (euler forward) 
     94         l_1st_euler = .TRUE.                 ! Set a Euler forward 1sr time-step 
    9595         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    9696         !                                    ! Initialization of ocean to zero 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/restart.F90

    r9839 r9863  
    88   !!            2.0  !  2006-07  (S. Masson)  use IOM for restart 
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
    10    !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D) 
    11    !!            3.7  !  2014-01  (G. Madec) suppression of curl and hdiv from the restart 
    12    !!             -   !  2014-12  (G. Madec) remove KPP scheme 
    13    !!---------------------------------------------------------------------- 
    14  
    15    !!---------------------------------------------------------------------- 
    16    !!   rst_opn    : open the ocean restart file 
    17    !!   rst_write  : write the ocean restart file 
    18    !!   rst_read   : read the ocean restart file 
    19    !!---------------------------------------------------------------------- 
    20    USE oce             ! ocean dynamics and tracers  
    21    USE dom_oce         ! ocean space and time domain 
    22    USE sbc_ice         ! only lk_si3  
    23    USE phycst          ! physical constants 
    24    USE eosbn2          ! equation of state            (eos bn2 routine) 
    25    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
     10   !!            - -  !  2010-10  (C. Ethe, G. Madec)  TRC-TRA merge (T-S in 4D) 
     11   !!            3.7  !  2014-01  (G. Madec)  suppression of curl and hdiv from the restart 
     12   !!             -   !  2014-12  (G. Madec)  remove KPP scheme 
     13   !!            4.0  !  2018-06  (G. Madec)  introduce l_1st_euler 
     14   !!---------------------------------------------------------------------- 
     15 
     16   !!---------------------------------------------------------------------- 
     17   !!   rst_opn       : open the ocean restart file in write mode 
     18   !!   rst_write     : write the ocean restart file 
     19   !!   rst_read_open : open the ocean restart file in read mode 
     20   !!   rst_read      : read the ocean restart file 
     21   !!---------------------------------------------------------------------- 
     22   USE oce            ! ocean dynamics and tracers  
     23   USE dom_oce        ! ocean space and time domain 
     24   USE sbc_ice        ! only lk_si3  
     25   USE phycst         ! physical constants 
     26   USE eosbn2         ! equation of state            (eos bn2 routine) 
     27   USE trdmxl_oce     ! ocean active mixed layer tracers trends variables 
     28   USE diurnal_bulk   !  
    2629   ! 
    27    USE in_out_manager  ! I/O manager 
    28    USE iom             ! I/O module 
    29    USE diurnal_bulk 
     30   USE in_out_manager ! I/O manager 
     31   USE iom            ! I/O module 
    3032 
    3133   IMPLICIT NONE 
    … …  
    3436   PUBLIC   rst_opn         ! routine called by step module 
    3537   PUBLIC   rst_write       ! routine called by step module 
     38   PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init 
    3639   PUBLIC   rst_read        ! routine called by istate module 
    37    PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init 
    3840 
    3941   !! * Substitutions 
    … …  
    144146      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
    145147      !!---------------------------------------------------------------------- 
    146                      IF(lwxios) CALL iom_swap(      cwxios_context          ) 
    147                      CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       , ldxios = lwxios)   ! dynamics time step 
    148  
    149       IF ( .NOT. ln_diurnal_only ) THEN 
    150                      CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub, ldxios = lwxios        )     ! before fields 
    151                      CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb, ldxios = lwxios        ) 
    152                      CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios ) 
    153                      CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios ) 
    154                      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb, ldxios = lwxios      ) 
    155                      ! 
    156                      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un, ldxios = lwxios        )     ! now fields 
    157                      CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn, ldxios = lwxios        ) 
    158                      CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios ) 
    159                      CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios ) 
    160                      CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn, ldxios = lwxios      ) 
    161                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      ) 
    162                   ! extra variable needed for the ice sheet coupling 
    163                   IF ( ln_iscpl ) THEN  
    164                      CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask, ldxios = lwxios ) ! need to extrapolate T/S 
    165                      CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask, ldxios = lwxios ) ! need to correct barotropic velocity 
    166                      CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask, ldxios = lwxios ) ! need to correct barotropic velocity 
    167                      CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask, ldxios = lwxios) ! need to correct barotropic velocity 
    168                      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )   ! need to compute temperature correction 
    169                      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation 
    170                      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation 
    171                      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n(:,:,:), ldxios = lwxios ) ! need to compute extrapolation if vvl 
    172                   END IF 
    173       ENDIF 
    174        
    175       IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios )   
    176       IF(lwxios) CALL iom_swap(      cxios_context          ) 
     148      IF( lwxios )   CALL iom_swap( cwxios_context ) 
     149          
     150      CALL    iom_rstput( kt, nitrst, numrow, 'rdt'    , rn_rdt           , ldxios = lwxios )   ! dynamics time step 
     151      ! 
     152      IF( .NOT. ln_diurnal_only ) THEN 
     153         CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub               , ldxios = lwxios )   ! before fields 
     154         CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb               , ldxios = lwxios ) 
     155         CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios ) 
     156         CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios ) 
     157         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb             , ldxios = lwxios ) 
     158         ! 
     159         CALL iom_rstput( kt, nitrst, numrow, 'un'     , un               , ldxios = lwxios )     ! now fields 
     160         CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn               , ldxios = lwxios ) 
     161         CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios ) 
     162         CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios ) 
     163         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn             , ldxios = lwxios ) 
     164         CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop             , ldxios = lwxios ) 
     165         ! 
     166         IF( ln_iscpl ) THEN          ! extra variable needed for the ice sheet coupling 
     167            CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask  , ldxios = lwxios )    ! need to extrapolate T/S 
     168            CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask  , ldxios = lwxios )    ! need to correct barotropic velocity 
     169            CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask  , ldxios = lwxios )    ! need to correct barotropic velocity 
     170            CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask , ldxios = lwxios )    ! need to correct barotropic velocity 
     171            CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n  , ldxios = lwxios )    ! need to compute temperature correction 
     172            CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n  , ldxios = lwxios )    ! need to compute bt conservation 
     173            CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n  , ldxios = lwxios )    ! need to compute bt conservation 
     174            CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n, ldxios = lwxios )    ! need to compute extrapolation if vvl 
     175         ENDIF 
     176      ENDIF 
     177      ! 
     178      IF( ln_diurnal )   CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios )   
     179      IF( lwxios     )   CALL iom_swap( cxios_context ) 
    177180      IF( kt == nitrst ) THEN 
    178          IF(.NOT.lwxios) THEN 
    179             CALL iom_close( numrow )     ! close the restart file (only at last time step) 
    180          ELSE 
    181             CALL iom_context_finalize(      cwxios_context          ) 
     181         IF( lwxios ) THEN   ;   CALL iom_context_finalize( cwxios_context ) 
     182         ELSE                ;   CALL iom_close( numrow )     ! close the restart file (only at last time step) 
    182183         ENDIF 
    183184!!gm         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE. 
    184185!!gm  not sure what to do here   ===>>>  ask to Sebastian 
    185186         lrst_oce = .FALSE. 
    186             IF( ln_rst_list ) THEN 
    187                nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1)) 
    188                nitrst = nstocklist( nrst_lst ) 
    189             ENDIF 
     187         IF( ln_rst_list ) THEN 
     188            nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1)) 
     189            nitrst = nstocklist( nrst_lst ) 
     190         ENDIF 
    190191      ENDIF 
    191192      ! 
    … …  
    202203      !!                the file has already been opened 
    203204      !!---------------------------------------------------------------------- 
    204       INTEGER        ::   jlibalt = jprstlib 
    205       LOGICAL        ::   llok 
    206       CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file 
     205      INTEGER       ::   jlibalt = jprstlib 
     206      LOGICAL       ::   llok 
     207      CHARACTER(lc) ::   clpath   ! full path to ocean output restart file 
    207208      !!---------------------------------------------------------------------- 
    208209      ! 
    … …  
    238239         ENDIF  
    239240      ENDIF 
    240  
     241      ! 
    241242   END SUBROUTINE rst_read_open 
    242243 
    … …  
    254255      REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d 
    255256      !!---------------------------------------------------------------------- 
    256  
     257      ! 
    257258      CALL rst_read_open           ! open restart for reading (if not already opened) 
    258259 
    … …  
    260261      IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN 
    261262         CALL iom_get( numror, 'rdt', zrdt, ldxios = lrxios ) 
    262          IF( zrdt /= rdt )   neuler = 0 
     263         IF( zrdt /= rn_rdt ) THEN 
     264            IF(lwp) WRITE( numout,*) 
     265            IF(lwp) WRITE( numout,*) 'rst_read:  rdt not equal to the read one' 
     266            IF(lwp) WRITE( numout,*) 
     267            IF(lwp) WRITE( numout,*) '      ==>>>   forced euler first time-step' 
     268            l_1st_euler =  .TRUE. 
     269         ENDIF 
    263270      ENDIF 
    264271 
    … …  
    266273      IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst, ldxios = lrxios )  
    267274      IF ( ln_diurnal_only ) THEN  
    268          IF(lwp) WRITE( numout, * ) & 
    269          &   "rst_read:- ln_diurnal_only set, setting rhop=rau0"  
     275         IF(lwp) WRITE( numout,*) 'rst_read: ln_diurnal_only set, setting rhop=rau0' 
    270276         rhop = rau0 
    271277         CALL iom_get( numror, jpdom_autoglo, 'tn'     , w3d, ldxios = lrxios )  
    … …  
    274280      ENDIF   
    275281       
    276       IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN 
     282      IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN 
    277283         CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub, ldxios = lrxios                )   ! before fields 
    278284         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb, ldxios = lrxios                ) 
    … …  
    281287         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios              ) 
    282288      ELSE 
    283          neuler = 0 
     289         l_1st_euler =  .TRUE.      ! before field not found, forced euler 1st time-step 
    284290      ENDIF 
    285291      ! 
    … …  
    295301      ENDIF 
    296302      ! 
    297       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0) 
    298          tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values 
    299          ub   (:,:,:)   = un   (:,:,:) 
    300          vb   (:,:,:)   = vn   (:,:,:) 
    301          sshb (:,:)     = sshn (:,:) 
    302          ! 
    303          IF( .NOT.ln_linssh ) THEN 
    304             DO jk = 1, jpk 
    305                e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    306             END DO 
    307          ENDIF 
    308          ! 
     303      IF( l_1st_euler ) THEN              ! Euler restart 
     304         tsb (:,:,:,:) = tsn (:,:,:,:)          ! all before fields set to now values 
     305         ub  (:,:,:)   = un  (:,:,:) 
     306         vb  (:,:,:)   = vn  (:,:,:) 
     307         sshb(:,:)     = sshn(:,:) 
     308         IF( .NOT.ln_linssh )   e3t_b(:,:,:) = e3t_n(:,:,:) 
    309309      ENDIF 
    310310      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynnxt.F90

    r9598 r9863  
    6464CONTAINS 
    6565 
    66    SUBROUTINE dyn_nxt ( kt ) 
     66   SUBROUTINE dyn_nxt( kt ) 
    6767      !!---------------------------------------------------------------------- 
    6868      !!                  ***  ROUTINE dyn_nxt  *** 
    … …  
    9292      !!               un,vn   now horizontal velocity of next time-step 
    9393      !!---------------------------------------------------------------------- 
    94       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     94      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    9797      INTEGER  ::   ikt          ! local integers 
    98       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars 
    99       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
     98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef   ! local scalars 
     99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf          !   -      - 
    100100      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    101101      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
    … …  
    132132            ! so that asselin contribution is removed at the same time  
    133133            DO jk = 1, jpkm1 
    134                un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 
    135                vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     134               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) ) * umask(:,:,jk) 
     135               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) ) * vmask(:,:,jk) 
    136136            END DO   
    137137         ENDIF 
    … …  
    152152!!$   Do we need a call to bdy_vol here?? 
    153153      ! 
    154       IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
    155          z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
    156          IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt 
    157          ! 
    158          !                                  ! Kinetic energy and Conversion 
    159          IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt ) 
    160          ! 
    161          IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends 
    162             zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 
    163             zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 
    164             CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter 
     154      IF( l_trddyn ) THEN             !* prepare the atf trend computation + some diagnostics 
     155         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )    ! Kinetic energy and Conversion 
     156         IF( ln_dyn_trd ) THEN                                       ! 3D output: total momentum trends 
     157            IF( ln_dynadv_vec ) THEN  
     158               zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_2dt 
     159               zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_2dt 
     160            ELSE 
     161               zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt 
     162               zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt 
     163            ENDIF 
     164            CALL iom_put( "utrd_tot", zua )                          ! total momentum trends, except the asselin filter 
    165165            CALL iom_put( "vtrd_tot", zva ) 
    166166         ENDIF 
    167          ! 
    168          zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter 
    169          zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the 
    170          !                                  !  computation of the asselin filter trends) 
     167         zua(:,:,:) = un(:,:,:)                    ! save the now velocity before the asselin filter 
     168         zva(:,:,:) = vn(:,:,:)                    ! (caution: the Asselin filter trends computation will be shifted by 1 timestep) 
    171169      ENDIF 
    172170 
    173171      ! Time filter and swap of dynamics arrays 
    174       ! ------------------------------------------ 
    175       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap 
     172      ! --------------------------------------- 
     173      IF( l_1st_euler ) THEN        !==  Euler at 1st time-step  ==!   (swap only) 
    176174         DO jk = 1, jpkm1 
    177175            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua 
    178176            vn(:,:,jk) = va(:,:,jk) 
    179177         END DO 
    180          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    181 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
     178         IF( .NOT.ln_linssh ) THEN                          ! e3._n <-- e3._a 
    182179            DO jk = 1, jpkm1 
    183 !               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    184 !               e3u_b(:,:,jk) = e3u_n(:,:,jk) 
    185 !               e3v_b(:,:,jk) = e3v_n(:,:,jk) 
    186                ! 
    187180               e3t_n(:,:,jk) = e3t_a(:,:,jk) 
    188181               e3u_n(:,:,jk) = e3u_a(:,:,jk) 
    189182               e3v_n(:,:,jk) = e3v_a(:,:,jk) 
    190183            END DO 
    191 !!gm BUG end 
    192          ENDIF 
    193                                                             !  
    194           
    195       ELSE                                             !* Leap-Frog : Asselin filter and swap 
     184         ENDIF 
     185         ! 
     186      ELSE                          !==  Leap-Frog  ==!   (Asselin filter and swap) 
     187         ! 
    196188         !                                ! =============! 
    197189         IF( ln_linssh ) THEN             ! Fixed volume ! 
    … …  
    213205         ELSE                             ! Variable volume ! 
    214206            !                             ! ================! 
    215             ! Before scale factor at t-points 
    216             ! (used as a now filtered scale factor until the swap) 
    217             ! ---------------------------------------------------- 
     207            ! Before scale factor at t-points   (used as a now filtered scale factor until the swap) 
    218208            DO jk = 1, jpkm1 
    219209               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 
    220210            END DO 
    221211            ! Add volume filter correction: compatibility with tracer advection scheme 
    222             ! => time filter + conservation correction (only at the first level) 
     212            !    => time filter + conservation correction (only at the first level) 
    223213            zcoef = atfp * rdt * r1_rau0 
    224214 
    … …  
    232222                           IF( jk <=  nk_rnf(ji,jj)  ) THEN 
    233223                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 
    234                                       &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
     224                                  &              * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 
    235225                           ENDIF 
    236                         ENDDO 
    237                      ENDDO 
    238                   ENDDO 
     226                        END DO 
     227                     END DO 
     228                  END DO 
    239229               ELSE 
    240230                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
    241231               ENDIF 
    242             END IF 
     232            ENDIF 
    243233 
    244234            IF ( ln_isf ) THEN   ! if ice shelf melting 
    … …  
    253243                  END DO 
    254244               END DO 
    255             END IF 
     245            ENDIF 
    256246            ! 
    257247            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    … …  
    322312         ENDIF 
    323313         ! 
    324       ENDIF ! neuler =/0 
     314      ENDIF    ! end Leap-Frog time stepping 
    325315      ! 
    326316      ! Set "now" and "before" barotropic velocities for next time step: 
    327       ! JC: Would be more clever to swap variables than to make a full vertical 
    328       ! integration 
     317      ! JC: Would be more clever to swap variables than to make a full vertical integration 
    329318      ! 
    330319      ! 
    … …  
    360349      ENDIF 
    361350      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum 
    362          zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 
    363          zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 
     351         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_2dt 
     352         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_2dt 
    364353         CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 
    365354      ENDIF 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg.F90

    r9598 r9863  
    7474      ! 
    7575      INTEGER  ::   ji, jj, jk                   ! dummy loop indices 
    76       REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars 
     76      REAL(wp) ::   zg_2, zintp, zgrau0r, zld   ! local scalars 
    7777      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice 
    7878      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90

    r9598 r9863  
    11MODULE dynspg_ts 
    2  
    3    !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !  
    4  
    52   !!====================================================================== 
    63   !!                   ***  MODULE  dynspg_ts  *** 
    … …  
    3532   USE sbcisf          ! ice shelf variable (fwfisf) 
    3633   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    37    USE dynadv    , ONLY: ln_dynadv_vec 
     34   USE dynadv   , ONLY : ln_dynadv_vec 
    3835   USE dynvor          ! vortivity scheme indicators 
    3936   USE phycst          ! physical constants 
    … …  
    8582   REAL(wp) ::   r1_2  = 0.5_wp           ! 
    8683 
     84   REAL(wp) ::   r1_2dt_b, r2dt_bf               ! local scalars 
     85    
    8786   !! * Substitutions 
    8887#  include "vectopt_loop_substitute.h90" 
    … …  
    151150      INTEGER  ::   ikbu, iktu, noffset   ! local integers 
    152151      INTEGER  ::   ikbv, iktv            !   -      - 
    153       REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars 
    154152      REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      - 
    155153      REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      - 
    … …  
    182180!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 
    183181      !                                         ! reciprocal of baroclinic time step  
    184       IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    185       ELSE                                        ;   z2dt_bf = 2.0_wp * rdt 
    186       ENDIF 
    187       r1_2dt_b = 1.0_wp / z2dt_bf  
     182      IF( l_1st_euler ) THEN   ;   r2dt_bf =          rdt 
     183      ELSE                     ;   r2dt_bf = 2.0_wp * rdt 
     184      ENDIF 
     185      r1_2dt_b = 1.0_wp / r2dt_bf  
    188186      ! 
    189187      ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart  
    … …  
    194192      ENDIF 
    195193      ! 
    196       IF( kt == nit000 ) THEN                   !* initialisation 
     194      IF( kt == nit000 ) THEN                   !* initialisation 1st time-step 
    197195         ! 
    198196         IF(lwp) WRITE(numout,*) 
    … …  
    201199         IF(lwp) WRITE(numout,*) 
    202200         ! 
    203          IF( neuler == 0 )   ll_init=.TRUE. 
    204          ! 
    205          IF( ln_bt_fw .OR. neuler == 0 ) THEN 
     201         IF( l_1st_euler )   ll_init = .TRUE. 
     202         ! 
     203         IF( ln_bt_fw .OR. l_1st_euler ) THEN 
    206204            ll_fw_start =.TRUE. 
    207205            noffset     = 0 
    … …  
    212210         ! Set averaging weights and cycle length: 
    213211         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
     212         ! 
     213      ELSEIF( kt == nit000 + 1 ) THEN           !* initialisation 2nd time-step 
     214         ! 
     215         IF( .NOT.ln_bt_fw .AND. l_1st_euler ) THEN 
     216            ll_fw_start = .FALSE. 
     217            CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
     218         ENDIF 
    214219         ! 
    215220      ENDIF 
    … …  
    340345         END SELECT 
    341346      ENDIF 
    342       ! 
    343       ! If forward start at previous time step, and centered integration,  
    344       ! then update averaging weights: 
    345       IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 
    346          ll_fw_start=.FALSE. 
    347          CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 
    348       ENDIF 
    349                            
     347      !                           
    350348      ! ----------------------------------------------------------------------------- 
    351349      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    … …  
    12031201         zwx(:,:) = un_adv(:,:) 
    12041202         zwy(:,:) = vn_adv(:,:) 
    1205          IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 
     1203         IF( .NOT.l_1st_euler ) THEN 
    12061204            un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 
    12071205            vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 
    … …  
    13051303      !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 
    13061304      !!---------------------------------------------------------------------- 
    1307       LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true. 
    1308       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true. 
    1309       INTEGER, INTENT(inout) :: jpit      ! cycle length     
    1310       REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights 
    1311                                                          zwgt2    ! Secondary weights 
    1312        
     1305      LOGICAL                       , INTENT(in   ) ::   ll_av          ! temporal averaging=.true. 
     1306      LOGICAL                       , INTENT(in   ) ::   ll_fw          ! forward time splitting =.true. 
     1307      INTEGER                       , INTENT(inout) ::   jpit           ! cycle length     
     1308      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, zwgt2   ! Primary & Secondary weights 
     1309      ! 
    13131310      INTEGER ::  jic, jn, ji                      ! temporary integers 
    13141311      REAL(wp) :: za1, za2 
    13151312      !!---------------------------------------------------------------------- 
    1316  
     1313      ! 
    13171314      zwgt1(:) = 0._wp 
    13181315      zwgt2(:) = 0._wp 
    1319  
     1316      ! 
    13201317      ! Set time index when averaged value is requested 
    1321       IF (ll_fw) THEN  
    1322          jic = nn_baro 
    1323       ELSE 
    1324          jic = 2 * nn_baro 
    1325       ENDIF 
    1326  
    1327       ! Set primary weights: 
    1328       IF (ll_av) THEN 
    1329            ! Define simple boxcar window for primary weights  
    1330            ! (width = nn_baro, centered around jic)      
     1318      IF ( ll_fw ) THEN   ;   jic =     nn_baro 
     1319      ELSE                ;   jic = 2 * nn_baro 
     1320      ENDIF 
     1321 
     1322      !                 !==  Set primary weights  ==! 
     1323      ! 
     1324      IF (ll_av) THEN         !* Define simple boxcar window for primary weights  
     1325         !                    !  (width = nn_baro, centered around jic)      
    13311326         SELECT CASE ( nn_bt_flt ) 
    1332               CASE( 0 )  ! No averaging 
    1333                  zwgt1(jic) = 1._wp 
    1334                  jpit = jic 
    1335  
    1336               CASE( 1 )  ! Boxcar, width = nn_baro 
    1337                  DO jn = 1, 3*nn_baro 
    1338                     za1 = ABS(float(jn-jic))/float(nn_baro)  
    1339                     IF (za1 < 0.5_wp) THEN 
    1340                       zwgt1(jn) = 1._wp 
    1341                       jpit = jn 
    1342                     ENDIF 
    1343                  ENDDO 
    1344  
    1345               CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
    1346                  DO jn = 1, 3*nn_baro 
    1347                     za1 = ABS(float(jn-jic))/float(nn_baro)  
    1348                     IF (za1 < 1._wp) THEN 
    1349                       zwgt1(jn) = 1._wp 
    1350                       jpit = jn 
    1351                     ENDIF 
    1352                  ENDDO 
    1353               CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
     1327         CASE( 0 )  ! No averaging 
     1328            zwgt1(jic) = 1._wp 
     1329            jpit = jic 
     1330            ! 
     1331         CASE( 1 )  ! Boxcar, width = nn_baro 
     1332            DO jn = 1, 3*nn_baro 
     1333               za1 = ABS(float(jn-jic))/float(nn_baro)  
     1334               IF ( za1 < 0.5_wp ) THEN 
     1335                  zwgt1(jn) = 1._wp 
     1336                  jpit = jn 
     1337               ENDIF 
     1338            END DO 
     1339            ! 
     1340         CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
     1341            DO jn = 1, 3*nn_baro 
     1342               za1 = ABS(float(jn-jic))/float(nn_baro)  
     1343                  IF ( za1 < 1._wp ) THEN 
     1344                     zwgt1(jn) = 1._wp 
     1345                     jpit = jn 
     1346                  ENDIF 
     1347               END DO 
     1348         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 
    13541349         END SELECT 
    1355  
    1356       ELSE ! No time averaging 
     1350         ! 
     1351      ELSE                    !* No time averaging 
    13571352         zwgt1(jic) = 1._wp 
    13581353         jpit = jic 
    13591354      ENDIF 
    13601355     
    1361       ! Set secondary weights 
     1356      !                 !==  Set secondary weights  ==! 
     1357      ! 
    13621358      DO jn = 1, jpit 
    1363         DO ji = jn, jpit 
    1364              zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
    1365         END DO 
     1359         DO ji = jn, jpit 
     1360            zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 
     1361         END DO 
    13661362      END DO 
    13671363 
    1368       ! Normalize weigths: 
    1369       za1 = 1._wp / SUM(zwgt1(1:jpit)) 
    1370       za2 = 1._wp / SUM(zwgt2(1:jpit)) 
     1364      !                 !==  Normalize weights  ==! 
     1365      ! 
     1366      za1 = 1._wp / SUM( zwgt1(1:jpit) ) 
     1367      za2 = 1._wp / SUM( zwgt2(1:jpit) ) 
    13711368      DO jn = 1, jpit 
    1372         zwgt1(jn) = zwgt1(jn) * za1 
    1373         zwgt2(jn) = zwgt2(jn) * za2 
     1369         zwgt1(jn) = zwgt1(jn) * za1 
     1370         zwgt2(jn) = zwgt2(jn) * za2 
    13741371      END DO 
    13751372      ! 
    … …  
    15391536      ! 
    15401537      !                             ! read restart when needed 
     1538!!gm what's happen when starting with an euler time-step BUT not from rest ? 
     1539!!   this case correspond to a restart with only now time-step available... 
    15411540      CALL ts_rst( nit000, 'READ' ) 
    15421541      ! 
    … …  
    15481547         CALL iom_set_rstw_var_active('vn_bf') 
    15491548         ! 
    1550          IF (.NOT.ln_bt_av) THEN 
     1549         IF ( .NOT.ln_bt_av ) THEN 
    15511550            CALL iom_set_rstw_var_active('sshbb_e') 
    15521551            CALL iom_set_rstw_var_active('ubb_e') 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90

    r9598 r9863  
    1111   !!---------------------------------------------------------------------- 
    1212   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing 
     13   !!       zdf_trd   : diagnose the zdf velocity trends and the KE dissipation trend  
     14!!gm                        ==>>> zdf_trd currently not used 
    1315   !!---------------------------------------------------------------------- 
    1416   USE oce            ! ocean dynamics and tracers variables 
    … …  
    2628   USE in_out_manager ! I/O manager 
    2729   USE lib_mpp        ! MPP library 
     30   USE iom             ! IOM library 
    2831   USE prtctl         ! Print control 
    2932   USE timing         ! Timing 
    … …  
    6770      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    6871      ! 
    69       INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    70       INTEGER  ::   iku, ikv           ! local integers 
    71       REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
    72       REAL(wp) ::   zzws, ze3va        !   -      - 
    73       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace  
    74       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      - 
     72      INTEGER  ::   ji, jj, jk            ! dummy loop indices 
     73      INTEGER  ::   iku, ikv              ! local integers 
     74      REAL(wp) ::   zzwi, ze3ua, z2dt_2   ! local scalars 
     75      REAL(wp) ::   zzws, ze3va           !   -      - 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace  
     77      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      - 
    7578      !!--------------------------------------------------------------------- 
    7679      ! 
    … …  
    8689         ENDIF 
    8790      ENDIF 
    88       !                             !* set time step 
    89       IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping) 
    90       ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog) 
    91       ENDIF 
     91      ! 
     92      z2dt_2 = r2dt * 0.5_wp        !* =rdt except in 1st Euler time step where it is equal to rdt/2 
     93      ! 
    9294      ! 
    9395      !                             !* explicit top/bottom drag case 
    … …  
    111113      ELSE                                      ! applied on thickness weighted velocity 
    112114         DO jk = 1, jpkm1 
    113             ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
    114                &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
    115             va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
    116                &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
     115            ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + r2dt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 
     116            va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + r2dt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
    117117         END DO 
    118118      ENDIF 
    … …  
    133133               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    134134               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    135                ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    136                va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     135               ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     136               va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
    137137            END DO 
    138138         END DO 
    … …  
    144144                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
    145145                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
    146                   ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    147                   va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     146                  ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     147                  va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    148148               END DO 
    149149            END DO 
    … …  
    153153      !              !==  Vertical diffusion on u  ==! 
    154154      ! 
    155       !                    !* Matrix construction 
    156       zdt = r2dt * 0.5 
    157       SELECT CASE( nldf_dyn ) 
    158       CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     155      SELECT CASE( nldf_dyn )    !* Matrix construction 
     156      ! 
     157      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu) 
    159158         DO jk = 1, jpkm1 
    160159            DO jj = 2, jpjm1  
    161160               DO ji = fs_2, fs_jpim1   ! vector opt. 
    162161                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
    163                   zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk ) + akzu(ji,jj,jk  ) )   & 
    164                      &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    165                   zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
    166                      &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     162                  zzwi = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) + akzu(ji,jj,jk  ) )   & 
     163                     &            / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     164                  zzws = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) )   & 
     165                     &            / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    167166                  zwi(ji,jj,jk) = zzwi 
    168167                  zws(ji,jj,jk) = zzws 
    … …  
    171170            END DO 
    172171         END DO 
    173       CASE DEFAULT               ! iso-level lateral mixing 
     172      CASE DEFAULT                  ! iso-level lateral mixing 
    174173         DO jk = 1, jpkm1 
    175174            DO jj = 2, jpjm1  
    176175               DO ji = fs_2, fs_jpim1   ! vector opt. 
    177176                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
    178                   zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
    179                   zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     177                  zzwi = - z2dt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     178                  zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
    180179                  zwi(ji,jj,jk) = zzwi 
    181180                  zws(ji,jj,jk) = zzws 
    … …  
    186185      END SELECT 
    187186      ! 
    188       DO jj = 2, jpjm1     !* Surface boundary conditions 
    189          DO ji = fs_2, fs_jpim1   ! vector opt. 
     187      DO jj = 2, jpjm1           !* Surface boundary conditions 
     188         DO ji = fs_2, fs_jpim1     ! vector opt. 
    190189            zwi(ji,jj,1) = 0._wp 
    191190            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    … …  
    204203               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
    205204               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
    206                zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     205               zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
    207206            END DO 
    208207         END DO 
    … …  
    213212                  iku = miku(ji,jj)       ! ocean top level at u- and v-points  
    214213                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
    215                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     214                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
    216215               END DO 
    217216            END DO 
    … …  
    245244         DO ji = fs_2, fs_jpim1   ! vector opt. 
    246245            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
    247             ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    248                &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
     246            ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rau0 ) * umask(ji,jj,1) 
    249247         END DO 
    250248      END DO 
    … …  
    272270      !              !==  Vertical diffusion on v  ==! 
    273271      ! 
    274       !                       !* Matrix construction 
    275       zdt = r2dt * 0.5 
    276       SELECT CASE( nldf_dyn ) 
    277       CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu) 
     272      !                       
     273      SELECT CASE( nldf_dyn )    !* Matrix construction 
     274      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu) 
    278275         DO jk = 1, jpkm1 
    279276            DO jj = 2, jpjm1    
    280277               DO ji = fs_2, fs_jpim1   ! vector opt. 
    281278                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    282                   zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk ) + akzv(ji,jj,jk  ) )   & 
    283                      &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    284                   zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
    285                      &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     279                  zzwi = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) + akzv(ji,jj,jk  ) )   & 
     280                     &            / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     281                  zzws = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) )   & 
     282                     &            / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    286283                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    287284                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
    … …  
    290287            END DO 
    291288         END DO 
    292       CASE DEFAULT               ! iso-level lateral mixing 
     289      CASE DEFAULT                  ! iso-level lateral mixing 
    293290         DO jk = 1, jpkm1 
    294291            DO jj = 2, jpjm1    
    295292               DO ji = fs_2, fs_jpim1   ! vector opt. 
    296293                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    297                   zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
    298                   zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     294                  zzwi = - z2dt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     295                  zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
    299296                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    300297                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
    … …  
    305302      END SELECT 
    306303      ! 
    307       DO jj = 2, jpjm1        !* Surface boundary conditions 
    308          DO ji = fs_2, fs_jpim1   ! vector opt. 
     304      DO jj = 2, jpjm1           !* Surface boundary conditions 
     305         DO ji = fs_2, fs_jpim1     ! vector opt. 
    309306            zwi(ji,jj,1) = 0._wp 
    310307            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    … …  
    322319               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    323320               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
    324                zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     321               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
    325322            END DO 
    326323         END DO 
    … …  
    330327                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    331328                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
    332                   zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
     329                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
    333330               END DO 
    334331            END DO 
    … …  
    362359         DO ji = fs_2, fs_jpim1   ! vector opt.           
    363360            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
    364             va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    365                &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
     361            va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rau0 ) * vmask(ji,jj,1) 
    366362         END DO 
    367363      END DO 
    … …  
    387383      END DO 
    388384      ! 
    389       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    390          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    391          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
     385      IF( l_trddyn ) THEN                        ! save the vertical diffusive trends for further diagnostics 
     386         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
     387            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_2dt - ztrdu(:,:,:) 
     388            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_2dt - ztrdv(:,:,:) 
     389         ELSE                                      ! applied on thickness weighted velocity 
     390            ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ztrdu(:,:,:) 
     391            ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ztrdv(:,:,:) 
     392         ENDIF 
    392393         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    393394         DEALLOCATE( ztrdu, ztrdv )  
    … …  
    401402   END SUBROUTINE dyn_zdf 
    402403 
     404!!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing 
     405 
     406   SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt ) 
     407      !!---------------------------------------------------------------------- 
     408      !!                  ***  ROUTINE zdf_trd  *** 
     409      !! 
     410      !! ** Purpose :   compute the trend due to the vert. momentum diffusion 
     411      !!              together with the Leap-Frog time stepping using an  
     412      !!              implicit scheme. 
     413      !! 
     414      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing 
     415      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf. 
     416      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise 
     417      !!               - update the after velocity with the implicit vertical mixing. 
     418      !!      This requires to solver the following system:  
     419      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 
     420      !!      with the following surface/top/bottom boundary condition: 
     421      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
     422      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
     423      !! 
     424      !! ** Action :   (ua,va)   after velocity  
     425      !!--------------------------------------------------------------------- 
     426      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptrdu, ptrdv   ! 3D work arrays use for zdf trends diag 
     427      INTEGER , INTENT(in   )                         ::   kt             ! ocean time-step index 
     428      ! 
     429      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     430      REAL(wp) ::   zzz              ! local scalar 
     431      REAL(wp) ::   zavmu, zavmum1   !   -      - 
     432      REAL(wp) ::   zavmv, zavmvm1   !   -      - 
     433      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2d    !  -      - 
     434      !!--------------------------------------------------------------------- 
     435      ! 
     436      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )   ! apply lateral boundary condition on (ua,va) 
     437      ! 
     438      ! 
     439      !                 !==  momentum trend due to vertical diffusion  ==! 
     440      ! 
     441      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
     442         ptrdu(:,:,:) = (              ua(:,:,:) -              ub(:,:,:) )                * r1_2dt - ptrdu(:,:,:) 
     443         ptrdv(:,:,:) = (              va(:,:,:) -              vb(:,:,:) )                * r1_2dt - ptrdv(:,:,:) 
     444      ELSE                                      ! applied on thickness weighted velocity 
     445         ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ptrdu(:,:,:) 
     446         ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ptrdv(:,:,:) 
     447      ENDIF 
     448      CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt ) 
     449      ! 
     450      ! 
     451      !                 !==  KE dissipation trend due to vertical diffusion  ==! 
     452      ! 
     453      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area 
     454         !                                ! due to v friction (v=vertical)  
     455         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points) 
     456         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of  
     457         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured). 
     458         !                                ! Note also that now e3 scale factors are used as after one are not computed ! 
     459         ! 
     460         CALL wrk_alloc(jpi,jpj,   z2d ) 
     461         z2d(:,:) = 0._wp 
     462         DO jk = 1, jpkm1 
     463            DO jj = 2, jpjm1 
     464               DO ji = 2, jpim1 
     465                  zavmu   = 0.5 * ( avm(ji+1,jj,jk) + avm(ji  ,jj,jk) ) 
     466                  zavmum1 = 0.5 * ( avm(ji  ,jj,jk) + avm(ji-1,jj,jk) ) 
     467                  zavmv   = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj  ,jk) ) 
     468                  zavmvm1 = 0.5 * ( avm(ji,jj  ,jk) + avm(ji,jj-1,jk) ) 
     469                
     470                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  & 
     471                     &   zavmu   * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / e3uw_n(ji  ,jj,jk) * wumask(ji  ,jj,jk)   & 
     472                     & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk)   & 
     473                     & + zavmv   * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / e3vw_n(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   & 
     474                     & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   & 
     475                     &                        ) 
     476!!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step  ===>>> use it ? 
     477!!                                                                                     No since in zdfshé only kz tke (or gls) is used 
     478!! 
     479!!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of 
     480!!gm     now by before shears, i.e. the source term of TKE (local positivity is not ensured). 
     481!!       CAUTION: requires to compute e3uw_a and e3vw_a !!! 
     482!                  z2d(ji,jj) = z2d(ji,jj)  + (                                                                                 & 
     483!                     &    avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) / e3uw_n(ji  ,jj,jk)                        & 
     484!                     &                     * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) ) / e3uw_a(ji  ,jj,jk) * wumask(ji  ,jj,jk)   & 
     485!                     &  + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk)                        & 
     486!                     &                       ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk)   & 
     487!                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) / e3vw_n(ji,jj  ,jk)                        & 
     488!                     &                       ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) ) / e3vw_a(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   & 
     489!                     &  + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk)                        & 
     490!                     &                       ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   & 
     491!                     &                       ) 
     492!!gm end 
     493               END DO 
     494            END DO 
     495         END DO 
     496         zzz= - 0.5_wp* rau0           ! caution sign minus here 
     497         z2d(:,:) = zzz * z2d(:,:)  
     498         CALL lbc_lnk( z2d,'T', 1. ) 
     499         CALL iom_put( 'dispkevfo', z2d ) 
     500      ENDIF 
     501      ! 
     502   END SUBROUTINE zdf_trd 
     503    
     504!!gm end not used 
     505 
    403506   !!============================================================================== 
    404507END MODULE dynzdf 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/sshwzv.F90

    r9598 r9863  
    6868      INTEGER, INTENT(in) ::   kt   ! time step 
    6969      !  
    70       INTEGER  ::   jk            ! dummy loop indice 
    71       REAL(wp) ::   z2dt, zcoef   ! local scalars 
     70      INTEGER  ::   jk         ! dummy loop indice 
     71      REAL(wp) ::   z1_2rau0   ! local scalars 
    7272      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace 
    7373      !!---------------------------------------------------------------------- 
    … …  
    8181      ENDIF 
    8282      ! 
    83       z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog) 
    84       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    85       zcoef = 0.5_wp * r1_rau0 
     83      z1_2rau0 = 0.5_wp * r1_rau0 
    8684 
    8785      !                                           !------------------------------! 
    8886      !                                           !   After Sea Surface Height   ! 
    8987      !                                           !------------------------------! 
    90       IF(ln_wd_il) THEN 
    91          CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
    92       ENDIF 
     88 
     89      IF(ln_wd_il)   CALL wad_lmt( sshb, z1_2rau0 * (emp_b(:,:) + emp(:,:)), r2dt ) 
    9390 
    9491      CALL div_hor( kt )                               ! Horizontal divergence 
    … …  
    10299      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    103100      !  
    104       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     101      ssha(:,:) = (  sshb(:,:) - r2dt * ( z1_2rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    105102      ! 
    106103#if defined key_agrif 
    … …  
    143140      ! 
    144141      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    145       REAL(wp) ::   z1_2dt       ! local scalars 
    146142      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv 
    147143      !!---------------------------------------------------------------------- 
    … …  
    159155      !                                           !     Now Vertical Velocity    ! 
    160156      !                                           !------------------------------! 
    161       z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog) 
    162       IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt 
    163157      ! 
    164158      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
    … …  
    180174            ! computation of w 
    181175            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    & 
    182                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk) 
     176               &                         + r1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk) 
    183177         END DO 
    184178         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
    … …  
    188182            ! computation of w 
    189183            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 & 
    190                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk) 
     184               &                         + r1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk) 
    191185         END DO 
    192186      ENDIF 
    … …  
    243237         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    244238      ENDIF 
    245       !              !==  Euler time-stepping: no filter, just swap  ==! 
    246       IF ( neuler == 0 .AND. kt == nit000 ) THEN 
    247          sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now) 
    248          ! 
    249       ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==! 
    250          !                                                  ! before <-- now filtered 
     239      !               
     240      IF ( l_1st_euler ) THEN    !==  Euler time-stepping  ==!   no filter, just swap 
     241         ! 
     242         sshn(:,:) = ssha(:,:)               ! now    <-- after  (before already = now) 
     243         ! 
     244      ELSE                       !==  Leap-Frog time-stepping  ==!   Asselin filter + swap 
     245         ! 
     246         !                                   ! before <-- now filtered 
    251247         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 
    252          IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed 
     248         IF( .NOT.ln_linssh ) THEN           ! before <-- with forcing removed 
    253249            zcoef = atfp * rdt * r1_rau0 
    254250            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
    … …  
    256252               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
    257253         ENDIF 
    258          sshn(:,:) = ssha(:,:)                              ! now <-- after 
     254         sshn(:,:) = ssha(:,:)               ! now <-- after 
    259255      ENDIF 
    260256      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/wet_dry.F90

    r9168 r9863  
    117117 
    118118 
    119    SUBROUTINE wad_lmt( sshb1, sshemp, z2dt ) 
     119   SUBROUTINE wad_lmt( sshb1, sshemp, p2dt ) 
    120120      !!---------------------------------------------------------------------- 
    121121      !!                  ***  ROUTINE wad_lmt  *** 
    … …  
    129129      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1        !!gm DOCTOR names: should start with p ! 
    130130      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sshemp 
    131       REAL(wp)                , INTENT(in   ) ::   z2dt 
     131      REAL(wp)                , INTENT(in   ) ::   p2dt 
    132132      ! 
    133133      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices 
    … …  
    220220                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,  jj-1) , 0._wp)  
    221221               ! 
    222                zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    223                zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
     222               zdep1 = (zzflxp + zzflxn) * p2dt / ztmp 
     223               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - p2dt * sshemp(ji,jj) 
    224224               ! 
    225225               IF( zdep1 > zdep2 ) THEN 
    226226                  wdmask(ji, jj) = 0._wp 
    227                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    228                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
     227                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * p2dt ) / ( zflxp(ji,jj) * p2dt ) 
     228                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * p2dt ) / ( zzflxp * p2dt ) 
    229229                  ! flag if the limiter has been used but stop flagging if the only 
    230230                  ! changes have zeroed the coefficient since further iterations will 
    … …  
    270270 
    271271 
    272    SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt ) 
     272   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, pdt ) 
    273273      !!---------------------------------------------------------------------- 
    274274      !!                  ***  ROUTINE wad_lmt  *** 
    … …  
    280280      !! ** Action  : - calculate flux limiter and W/D flag 
    281281      !!---------------------------------------------------------------------- 
    282       REAL(wp)                , INTENT(in   ) ::   rdtbt    ! ocean time-step index 
     282      REAL(wp)                , INTENT(in   ) ::   pdt    ! external mode time-step [s] 
    283283      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc   
    284284      ! 
    285285      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices 
    286286      INTEGER  ::   jflag               ! local integer 
    287       REAL(wp) ::   z2dt 
    288287      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars 
    289288      REAL(wp) ::   zzflxp, zzflxn      ! local scalars 
    … …  
    298297      jflag  = 0 
    299298      zdepwd = 50._wp   ! maximum depth that ocean cells can have W/D processes 
    300       ! 
    301       z2dt = rdtbt    
    302299      ! 
    303300      zflxp(:,:)   = 0._wp 
    … …  
    347344                  &   + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
    348345           
    349                zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    350                zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
     346               zdep1 = (zzflxp + zzflxn) * pdt / ztmp 
     347               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - pdt * zssh_frc(ji,jj) 
    351348           
    352349               IF(zdep1 > zdep2) THEN 
    353                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    354                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt ) 
     350                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * pdt ) / ( zflxp(ji,jj) * pdt ) 
     351                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * pdt ) / ( zzflxp * pdt ) 
    355352                 ! flag if the limiter has been used but stop flagging if the only 
    356353                 ! changes have zeroed the coefficient since further iterations will 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traadv.F90

    r9598 r9863  
    9292      IF( ln_timing )   CALL timing_start('tra_adv') 
    9393      ! 
    94       !                                          ! set time step 
    95       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =         rdt   ! at nit000             (Euler) 
    96       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp * rdt   ! at nit000 or nit000+1 (Leapfrog) 
    97       ENDIF 
    98       ! 
    9994      !                                         !==  effective transport  ==! 
    10095      zun(:,:,jpk) = 0._wp 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf.F90

    r9598 r9863  
    5555      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5656      !! 
    57       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
     57      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd   ! 4D workspace 
    5858      !!---------------------------------------------------------------------- 
    5959      ! 
    … …  
    6161      ! 
    6262      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    63          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )  
    64          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
    65          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     63         ALLOCATE( ztrd(jpi,jpj,jpk,jpts) )  
     64         ztrd(:,:,:,:) = tsa(:,:,:,:)  
    6665      ENDIF 
    6766      ! 
    … …  
    7877      ! 
    7978      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    80          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    81          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    82          CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    83          CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 
    84          DEALLOCATE( ztrdt, ztrds )  
     79         ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:) 
     80         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrd(:,:,:,jp_tem) ) 
     81         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrd(:,:,:,jp_sal) ) 
     82         DEALLOCATE( ztrd )  
    8583      ENDIF 
    8684      !                                        !* print mean trends (used for debugging) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf_iso.F90

    r9779 r9863  
    108108      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
    109109      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
    110       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
     110      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      - 
    111111      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d 
    112112      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw  
    … …  
    127127      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    128128         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    129       ! 
    130       !                                            ! set time step size (Euler/Leapfrog) 
    131       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    132       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
    133       ENDIF 
    134       z1_2dt = 1._wp / z2dt 
    135129      ! 
    136130      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    … …  
    191185                     DO ji = 1, fs_jpim1 
    192186                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
    193                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    194                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     187                        zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     188                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt 
    195189                     END DO 
    196190                  END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf_triad.F90

    r9598 r9863  
    8585      INTEGER  ::  ip,jp,kp         ! dummy loop indices 
    8686      INTEGER  ::  ierr            ! local integer 
    87       REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars 
    88       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      - 
    89       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      - 
     87      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
     88      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
     89      REAL(wp) ::  zcoef0, ze3w_2, zsign         !   -      - 
    9090      ! 
    9191      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv 
    … …  
    110110      l_hst = .FALSE. 
    111111      l_ptr = .FALSE. 
    112       IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE.  
    113       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    114          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    115       ! 
    116       !                                                        ! set time step size (Euler/Leapfrog) 
    117       IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    118       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
     112      IF( cdtype == 'TRA' ) THEN 
     113         IF ( ln_diaptr )                                                 l_ptr = .TRUE.  
     114         IF ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
     115            & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")      )   l_hst = .TRUE. 
    119116      ENDIF 
    120       z1_2dt = 1._wp / z2dt 
    121117      ! 
    122118      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    … …  
    202198                     DO ji = 1, fs_jpim1 
    203199                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 
    204                         zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    205                         akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     200                        zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     201                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt 
    206202                     END DO 
    207203                  END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/tranpc.F90

    r9598 r9863  
    6565      LOGICAL  ::   l_bottom_reached, l_column_treated 
    6666      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
    67       REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     67      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw 
    6868      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0) 
    6969      REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point... 
    … …  
    301301         ! 
    302302         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
    303             z1_r2dt = 1._wp / (2._wp * rdt) 
    304             ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt 
    305             ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt 
     303            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * r1_2dt 
     304            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * r1_2dt 
    306305            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
    307306            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/tranxt.F90

    r9598 r9863  
    111111      IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries 
    112112  
    113       ! set time step size (Euler/Leapfrog) 
    114       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler) 
    115       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog) 
    116       ENDIF 
    117  
    118       ! trends computation initialisation 
    119       IF( l_trdtra )   THEN                     
     113      IF( l_trdtra )   THEN               ! trends computation initialisation 
    120114         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    121115         ztrdt(:,:,jpk) = 0._wp 
    … …  
    142136      ENDIF 
    143137 
    144       IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap) 
     138      IF( l_1st_euler ) THEN        ! Euler time-stepping at first time-step (only swap) 
    145139         DO jn = 1, jpts 
    146140            DO jk = 1, jpkm1 
    … …  
    156150         END IF 
    157151         ! 
    158       ELSE                                            ! Leap-Frog + Asselin filter time stepping 
     152      ELSE                          ! Leap-Frog + Asselin filter time stepping 
    159153         ! 
    160154         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface  
    … …  
    163157         ENDIF 
    164158         ! 
    165          CALL lbc_lnk_multi( tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., & 
    166                   &          tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., & 
    167                   &          tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  ) 
     159!!gm 
     160!         CALL lbc_lnk_multi( tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., & 
     161!            &                tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., & 
     162!            &                tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  ) 
     163!!gm 
     164         CALL lbc_lnk_multi( tsb, 'T', 1., tsn, 'T', 1., tsa, 'T', 1.  ) 
     165!!gm 
    168166         ! 
    169167      ENDIF      

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traqsr.F90

    r9598 r9863  
    133133      !                         !-----------------------------------! 
    134134      IF( kt == nit000 ) THEN          !==  1st time step  ==! 
    135 !!gm case neuler  not taken into account.... 
    136          IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart 
     135         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart 
    137136            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file' 
    138137            z1_2 = 0.5_wp 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/trazdf.F90

    r9598 r9863  
    5252      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    5353      ! 
    54       INTEGER  ::   jk   ! Dummy loop indices 
    55       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace 
     54      INTEGER  ::   jk, jts   ! Dummy loop indices 
     55      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   ztrd   ! 4D workspace 
    5656      !!--------------------------------------------------------------------- 
    5757      ! 
    … …  
    6464      ENDIF 
    6565      ! 
    66       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping) 
    67       ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog) 
    68       ENDIF 
    69       ! 
    7066      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
    71          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    72          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    73          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     67         ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 
     68         ztrd(:,:,:,:) = tsa(:,:,:,:) 
    7469      ENDIF 
    7570      ! 
    … …  
    8580 
    8681      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    87          DO jk = 1, jpkm1 
    88             ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 
    89                &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 
    90             ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 
    91               &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 
     82         DO jts = 1, jpts 
     83            DO jk = 1, jpkm1 
     84               ztrd(:,:,jk,jts) = ( ( tsa(:,:,jk,jts)*e3t_a(:,:,jk) - tsb(:,:,jk,jts)*e3t_b(:,:,jk) ) / (e3t_n(:,:,jk)*r2dt) )   & 
     85                  &            - ztrd(:,:,jk,jts) 
     86            END DO 
    9287         END DO 
    9388!!gm this should be moved in trdtra.F90 and done on all trends 
    94          CALL lbc_lnk_multi( ztrdt, 'T', 1. , ztrds, 'T', 1. ) 
     89         CALL lbc_lnk( ztrd, 'T', 1. ) 
    9590!!gm 
    96          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    97          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    98          DEALLOCATE( ztrdt , ztrds ) 
     91         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrd(:,:,:,jp_tem) ) 
     92         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrd(:,:,:,jp_sal) ) 
     93         DEALLOCATE( ztrd ) 
    9994      ENDIF 
    10095      !                                          ! print mean trends (used for debugging) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRD/trdtra.F90

    r9598 r9863  
    237237      INTEGER                   , INTENT(in   ) ::   kt      ! time step 
    238238      !!---------------------------------------------------------------------- 
    239  
    240       IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdt (restart with Euler time stepping) 
    241       ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdt (leapfrog) 
    242       ENDIF 
    243239 
    244240      !                   ! 3D output of tracers trends using IOM interface 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/nemogcm.F90

    r9780 r9863  
    151151      !                            !==   time stepping   ==! 
    152152      !                            !-----------------------! 
     153      ! 
     154      !                                               !== set the model time-step  ==! 
     155      ! 
     156      IF( l_1st_euler ) THEN   ;   r2dt =         rn_rdt   ;   l_1st_euler = .TRUE.    ! start or restart with Euler 1st time-step 
     157      ELSE                     ;   r2dt = 2._wp * rn_rdt   ;   l_1st_euler = .FALSE.   ! restart with leapfrog 
     158      ENDIF 
     159      r1_2dt = 1._wp / r2dt 
     160      ! NB: if l_1st_euler=T, r2dt will be set to 2*rdt at the end of the 1st time-step (in step.F90) 
     161      !     Done here (not in domain.F90) as in ASM initialization an Euler 1st time step can be forced 
     162      ! 
     163      ! 
    153164      istp = nit000 
    154165      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/step.F90

    r9780 r9863  
    3434 
    3535   !!---------------------------------------------------------------------- 
    36    !!   stp             : OPA system time-stepping 
    37    !!---------------------------------------------------------------------- 
    38    USE step_oce         ! time stepping definition modules 
     36   !!   stp           : OPA system time-stepping 
     37   !!---------------------------------------------------------------------- 
     38   USE step_oce       ! time stepping definition modules 
    3939   ! 
    40    USE iom              ! xIOs server 
     40   USE iom            ! xIOs server 
    4141 
    4242   IMPLICIT NONE 
    … …  
    323323#endif 
    324324      ! 
     325      IF( l_1st_euler ) THEN 
     326         r2dt   = 2._wp * rn_rdt                            ! recover Leap-frog time-step 
     327         r1_2dt = 1._wp / r2dt 
     328         l_1st_euler = .FALSE. 
     329      ENDIF 
     330      ! 
    325331      IF( ln_timing )   CALL timing_stop('stp') 
    326332      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OFF/nemogcm.F90

    r9751 r9863  
    100100      !                            !==   time stepping   ==! 
    101101      !                            !-----------------------! 
     102      !                                               !== set the model time-step  ==! 
     103      ! 
     104      IF( l_1st_euler ) THEN   ;   r2dt =         rn_rdt   ;   l_1st_euler = .TRUE.    ! start or restart with Euler 1st time-step 
     105      ELSE                     ;   r2dt = 2._wp * rn_rdt   ;   l_1st_euler = .FALSE.   ! restart with leapfrog 
     106      ENDIF 
     107      r1_2dt = 1._wp / r2dt 
     108      ! NB: if l_1st_euler=T, r2dt will be set to 2*rdt at the end of the 1st time-step (in step.F90) 
     109      !     Done here (not in domain.F90) as in ASM initialization an Euler 1st time step can be forced 
     110      ! 
     111      ! 
    102112      istp = nit000 
    103113      ! 
    … …  
    115125                                CALL trc_stp    ( istp )         ! time-stepping 
    116126                                CALL stp_ctl    ( istp, indic )  ! Time loop: control and print 
     127         IF( l_1st_euler ) THEN 
     128            r2dt   = 2._wp * rn_rdt                              ! recover Leap-frog time-step 
     129            r1_2dt = 1._wp / r2dt 
     130            l_1st_euler = .FALSE. 
     131         ENDIF 
     132         ! 
    117133         istp = istp + 1 
    118134      END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/PISCES/P2Z/p2zexp.F90

    r9788 r9863  
    44   !! TOP :   LOBSTER Compute loss of organic matter in the sediments 
    55   !!====================================================================== 
    6    !! History :    -   !  1999    (O. Aumont, C. Le Quere)  original code 
    7    !!              -   !  2001-05 (O. Aumont, E. Kestenare) add sediment computations 
    8    !!             1.0  !  2005-06 (A.-S. Kremeur) new temporal integration for sedpoc 
    9    !!             2.0  !  2007-12  (C. Deltel, G. Madec)  F90 
    10    !!             3.5  !  2012-03  (C. Ethe)  Merge PISCES-LOBSTER 
    11    !!---------------------------------------------------------------------- 
    12    !!   p2z_exp        :  Compute loss of organic matter in the sediments 
    13    !!---------------------------------------------------------------------- 
    14    USE oce_trc         ! 
    15    USE trc 
    16    USE sms_pisces 
    17    USE p2zsed 
    18    USE lbclnk 
    19    USE prtctl_trc      ! Print control for debbuging 
    20    USE trd_oce 
    21    USE trdtrc 
    22    USE iom 
     6   !! History :   -   !  1999    (O. Aumont, C. Le Quere)  original code 
     7   !!             -   !  2001-05 (O. Aumont, E. Kestenare) add sediment computations 
     8   !!            1.0  !  2005-06 (A.-S. Kremeur) new temporal integration for sedpoc 
     9   !!            2.0  !  2007-12  (C. Deltel, G. Madec)  F90 
     10   !!            3.5  !  2012-03  (C. Ethe)  Merge PISCES-LOBSTER 
     11   !!---------------------------------------------------------------------- 
     12   !!   p2z_exp       :  Compute loss of organic matter in the sediments 
     13   !!---------------------------------------------------------------------- 
     14   USE oce_trc        ! 
     15   USE trc            ! 
     16   USE sms_pisces     ! 
     17   USE p2zsed         ! 
     18   USE lbclnk         ! 
     19   USE prtctl_trc     ! Print control for debbuging 
     20   USE trd_oce        ! 
     21   USE trdtrc         ! 
     22   USE iom            ! 
    2323 
    2424   IMPLICIT NONE 
    … …  
    3030 
    3131   ! 
    32    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   dminl     !: fraction of sinking POC released in sediments 
    33    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   dmin3     !: fraction of sinking POC released at each level 
    34    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   sedpocb   !: mass of POC in sediments 
    35    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   sedpocn   !: mass of POC in sediments 
    36    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   cmask     !: Coastal mask area 
    37    REAL(wp)                                ::   areacot   !: surface coastal area 
     32   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   dminl     ! fraction of sinking POC released in sediments 
     33   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   dmin3     ! fraction of sinking POC released at each level 
     34   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   sedpocb   ! mass of POC in sediments 
     35   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   sedpocn   ! mass of POC in sediments 
     36   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   cmask     ! Coastal mask area 
     37   REAL(wp)                                ::   areacot   ! surface coastal area 
    3838 
    3939   !! * Substitutions 
    … …  
    5959      !!              COLUMN BELOW THE SURFACE LAYER. 
    6060      !!--------------------------------------------------------------------- 
    61       !! 
    62       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
     61      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index       
    6362      !! 
    6463      INTEGER  ::   ji, jj, jk, jl, ikt 
    6564      REAL(wp) ::   zgeolpoc, zfact, zwork, ze3t, zsedpocd, zmaskt 
    6665      REAL(wp), DIMENSION(jpi,jpj)   ::  zsedpoca 
    67       CHARACTER (len=25) :: charout 
     66      CHARACTER (len=25) ::   charout 
    6867      !!--------------------------------------------------------------------- 
    6968      ! 
    … …  
    7271      IF( kt == nittrc000 )   CALL p2z_exp_init 
    7372 
    74       zsedpoca(:,:) = 0. 
    75  
    76  
    77       ! VERTICAL DISTRIBUTION OF NEWLY PRODUCED BIOGENIC 
    78       ! POC IN THE WATER COLUMN 
     73      zsedpoca(:,:) = 0._wp 
     74 
     75 
     76      ! VERTICAL DISTRIBUTION OF NEWLY PRODUCED BIOGENIC POC IN THE WATER COLUMN 
    7977      ! (PARTS OF NEWLY FORMED MATTER REMAINING IN THE DIFFERENT 
    8078      ! LAYERS IS DETERMINED BY DMIN3 DEFINED IN sms_p2z.F90 
    … …  
    9391    
    9492 
    95       zgeolpoc = 0.e0         !     Initialization 
    96       ! Release of nutrients from the "simple" sediment 
    97       DO jj = 2, jpjm1 
     93      zgeolpoc = 0._wp        !     Initialization 
     94      DO jj = 2, jpjm1           ! Release of nutrients from the "simple" sediment 
    9895         DO ji = fs_2, fs_jpim1 
    9996            ikt = mbkt(ji,jj)  
    … …  
    121118      ! Time filter and swap of arrays 
    122119      ! ------------------------------ 
    123       IF( neuler == 0 .AND. kt == nittrc000 ) THEN        ! Euler time-stepping at first time-step 
    124         !                                             ! (only swap) 
     120      IF( l_1st_euler ) THEN        ! Euler time-stepping at first time-step   (only swap) 
    125121        sedpocn(:,:) = zsedpoca(:,:) 
    126122        !                                               
    127       ELSE 
     123      ELSE                          ! Leap-Frog + Asselin filter 
    128124        ! 
    129125        DO jj = 1, jpj 
    130126           DO ji = 1, jpi 
    131               zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj)      ! time laplacian on tracers 
     127              zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj)     ! time laplacian on tracers 
    132128              sedpocb(ji,jj) = sedpocn(ji,jj) + atfp * zsedpocd                     ! sedpocb <-- filtered sedpocn 
    133               sedpocn(ji,jj) = zsedpoca(ji,jj)                                       ! sedpocn <-- sedpoca 
     129              sedpocn(ji,jj) = zsedpoca(ji,jj)                                      ! sedpocn <-- sedpoca 
    134130           END DO 
    135131        END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/PISCES/P2Z/p2zsms.F90

    r9598 r9863  
    44   !! TOP :   Time loop of LOBSTER model 
    55   !!====================================================================== 
    6    !! History :   1.0  !            M. Levy 
    7    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  revised architecture 
     6   !! History :  1.0  !            M. Levy 
     7   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  revised architecture 
    88   !!---------------------------------------------------------------------- 
    99 
    … …  
    1111   !!   p2zsms        :  Time loop of passive tracers sms 
    1212   !!---------------------------------------------------------------------- 
    13    USE oce_trc          ! 
    14    USE trc 
    15    USE sms_pisces 
    16    USE p2zbio 
    17    USE p2zopt 
    18    USE p2zsed 
    19    USE p2zexp 
    20    USE trd_oce 
    21    USE trdtrc_oce 
    22    USE trdtrc 
    23    USE trdmxl_trc 
     13   USE oce_trc        ! 
     14   USE trc            ! 
     15   USE sms_pisces     ! 
     16   USE p2zbio         ! 
     17   USE p2zopt         ! 
     18   USE p2zsed         ! 
     19   USE p2zexp         ! 
     20   USE trd_oce        ! 
     21   USE trdtrc_oce     ! 
     22   USE trdtrc         ! 
     23   USE trdmxl_trc     ! 
    2424 
    2525   IMPLICIT NONE 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/PISCES/P4Z/p4zsms.F90

    r9751 r9863  
    44   !! TOP :   PISCES Source Minus Sink manager 
    55   !!====================================================================== 
    6    !! History :   1.0  !  2004-03 (O. Aumont) Original code 
    7    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
     6   !! History :  1.0  !  2004-03 (O. Aumont) Original code 
     7   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
    88   !!---------------------------------------------------------------------- 
    9    !!   p4z_sms        : Time loop of passive tracers sms 
     9 
    1010   !!---------------------------------------------------------------------- 
    11    USE oce_trc         ! shared variables between ocean and passive tracers 
    12    USE trc             ! passive tracers common variables  
    13    USE trcdta          !  
    14    USE sms_pisces      ! PISCES Source Minus Sink variables 
    15    USE p4zbio          ! Biological model 
    16    USE p4zche          ! Chemical model 
    17    USE p4zlys          ! Calcite saturation 
    18    USE p4zflx          ! Gas exchange 
    19    USE p4zsbc          ! External source of nutrients 
    20    USE p4zsed          ! Sedimentation 
    21    USE p4zint          ! time interpolation 
    22    USE p4zrem          ! remineralisation 
    23    USE iom             ! I/O manager 
    24    USE trd_oce         ! Ocean trends variables 
    25    USE trdtrc          ! TOP trends variables 
    26    USE sedmodel        ! Sediment model 
    27    USE prtctl_trc      ! print control for debugging 
     11   !!   p4z_sms       : Time loop of passive tracers sms 
     12   !!   p4z_sms_init  : initialisation 
     13   !!   p4z_rst       : Read or write variables in restart file 
     14   !!   p4z_dmp       : Relaxation of some tracers 
     15   !!   p4z_chk_mass  : mass conservation check 
     16   !!---------------------------------------------------------------------- 
     17   USE oce_trc        ! shared variables between ocean and passive tracers 
     18   USE trc            ! passive tracers common variables  
     19   USE trcdta         !  
     20   USE sms_pisces     ! PISCES Source Minus Sink variables 
     21   USE p4zbio         ! Biological model 
     22   USE p4zche         ! Chemical model 
     23   USE p4zlys         ! Calcite saturation 
     24   USE p4zflx         ! Gas exchange 
     25   USE p4zsbc         ! External source of nutrients 
     26   USE p4zsed         ! Sedimentation 
     27   USE p4zint         ! time interpolation 
     28   USE p4zrem         ! remineralisation 
     29   USE trd_oce        ! Ocean trends variables 
     30   USE trdtrc         ! TOP trends variables 
     31   USE sedmodel       ! Sediment model 
     32   ! 
     33   USE iom            ! I/O manager 
     34   USE prtctl_trc     ! print control for debugging 
    2835 
    2936   IMPLICIT NONE 
    … …  
    3744   REAL(wp) ::   xfact1, xfact2, xfact3 
    3845 
    39    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     ! Array used to indicate negative tracer values 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr   ! Array used to indicate negative tracer values 
    4047 
    4148   !!---------------------------------------------------------------------- 
    … …  
    95102      ENDIF 
    96103 
    97       IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN 
     104      IF( l_1st_euler .OR. ln_top_euler ) THEN 
    98105         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter 
    99106            trb(:,:,:,jn) = trn(:,:,:,jn) 
    … …  
    277284         IF(lwp) WRITE(numout,*) 
    278285         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model ' 
    279          IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
     286         IF(lwp) WRITE(numout,*) ' ~~~~~~~' 
    280287         !  
    281288         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN 
    … …  
    407414      !! 
    408415      !!--------------------------------------------------------------------- 
    409       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
     416      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     417      ! 
    410418      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot 
    411419      CHARACTER(LEN=100)   ::   cltxt 
    412       INTEGER :: jk 
    413       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwork 
     420      INTEGER ::   jk 
     421      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwork 
    414422      !!---------------------------------------------------------------------- 
    415423      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/TRP/trcnxt.F90

    r9598 r9863  
    2424   !!   'key_top'                                                TOP models 
    2525   !!---------------------------------------------------------------------- 
    26    !!   trc_nxt     : time stepping on passive tracers 
    27    !!---------------------------------------------------------------------- 
    28    USE oce_trc         ! ocean dynamics and tracers variables 
    29    USE trc             ! ocean passive tracers variables 
    30    USE trd_oce 
    31    USE trdtra 
    32    USE tranxt 
    33    USE bdy_oce   , ONLY: ln_bdy 
    34    USE trcbdy          ! BDY open boundaries 
     26 
     27   !!---------------------------------------------------------------------- 
     28   !!   trc_nxt       : time stepping on passive tracers 
     29   !!---------------------------------------------------------------------- 
     30   USE oce_trc        ! ocean dynamics and tracers variables 
     31   USE trc            ! ocean passive tracers variables 
     32   USE trd_oce        !  
     33   USE trdtra         ! 
     34   USE tranxt         ! 
     35   USE bdy_oce , ONLY : ln_bdy 
     36   USE trcbdy         ! BDY open boundaries 
    3537# if defined key_agrif 
    3638   USE agrif_top_interp 
    3739# endif 
    3840   ! 
    39    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    40    USE prtctl_trc      ! Print control for debbuging 
     41   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     42   USE prtctl_trc     ! Print control for debbuging 
    4143 
    4244   IMPLICIT NONE 
    … …  
    8183      ! 
    8284      INTEGER  ::   jk, jn   ! dummy loop indices 
    83       REAL(wp) ::   zfact            ! temporary scalar 
     85      REAL(wp) ::   zfact    ! local scalar 
    8486      CHARACTER (len=22) :: charout 
    8587      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace 
    … …  
    99101      CALL lbc_lnk( tra(:,:,:,:), 'T', 1. )    
    100102 
    101       IF( ln_bdy )  CALL trc_bdy( kt ) 
     103      IF( ln_bdy )   CALL trc_bdy( kt ) 
    102104 
    103105      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application 
    104106         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) ) 
    105          ztrdt(:,:,:,:)  = trn(:,:,:,:) 
     107         ztrdt(:,:,:,:) = trn(:,:,:,:) 
    106108      ENDIF 
    107109      !                                ! Leap-Frog + Asselin filter time stepping 
    108       IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap) 
     110      IF( l_1st_euler .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap) 
    109111         DO jn = 1, jptra 
    110112            DO jk = 1, jpkm1 
    … …  
    208210                  ztc_f  = ztc_n  + atfp * ztc_d 
    209211                  ! 
    210                   IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level  
     212                  IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! top ocean level  
    211213                     ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   )  
    212214                     ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) ) 
    213215                  ENDIF 
    214  
    215                   ze3t_f = 1.e0 / ze3t_f 
    216                   trb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     216                  ! 
     217                  trb(ji,jj,jk,jn) = ztc_f / ze3t_f       ! ptb <-- ptn filtered 
    217218                  trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn)     ! ptn <-- pta 
    218                   ! 
    219219               END DO 
    220220            END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/trc.F90

    r9598 r9863  
    6464   CHARACTER(len = 256), PUBLIC ::   cn_trcrst_outdir   !: restart output directory 
    6565   REAL(wp)            , PUBLIC ::   rdttrc             !: passive tracer time step 
    66    REAL(wp)            , PUBLIC ::   r2dttrc            !: = 2*rdttrc except at nit000 (=rdttrc) if neuler=0 
     66   REAL(wp)            , PUBLIC ::   r2dttrc            !: = 2*rdttrc except at nit000 (=rdttrc) if l_top_euler=T 
     67   REAL(wp)            , PUBLIC ::   r1_2dttrc          !: = 1/rdttrc 
    6768   LOGICAL             , PUBLIC ::   ln_top_euler       !: boolean term for euler integration  
     69   LOGICAL             , PUBLIC ::   l_top_euler        !: boolean term for euler integration  
    6870   LOGICAL             , PUBLIC ::   ln_trcdta          !: Read inputs data from files 
    6971   LOGICAL             , PUBLIC ::   ln_trcdmp          !: internal damping flag 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/trcrst.F90

    r9598 r9863  
    44   !! TOP :   Manage the passive tracer restart 
    55   !!====================================================================== 
    6    !! History :    -   !  1991-03  ()  original code 
    7    !!             1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90 
    8    !!              -   !  2005-10 (C. Ethe) print control 
    9    !!             2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture 
     6   !! History :   -   !  1991-03  ()  original code 
     7   !!            1.0  !  2005-03 (O. Aumont, A. El Moussaoui) F90 
     8   !!             -   !  2005-10 (C. Ethe) print control 
     9   !!            2.0  !  2005-10 (C. Ethe, G. Madec) revised architecture 
    1010   !!---------------------------------------------------------------------- 
    1111#if defined key_top 
    … …  
    1313   !!   'key_top'                                                TOP models 
    1414   !!---------------------------------------------------------------------- 
    15    !!---------------------------------------------------------------------- 
    16    !!   trc_rst        : Restart for passive tracer 
    17    !!   trc_rst_opn    : open  restart file 
    18    !!   trc_rst_read   : read  restart file 
    19    !!   trc_rst_wri    : write restart file 
    20    !!---------------------------------------------------------------------- 
    21    USE oce_trc 
    22    USE trc 
    23    USE iom 
    24    USE daymod 
     15    
     16   !!---------------------------------------------------------------------- 
     17   !!   trc_rst       : Restart for passive tracer 
     18   !!   trc_rst_opn   : open  restart file 
     19   !!   trc_rst_read  : read  restart file 
     20   !!   trc_rst_wri   : write restart file 
     21   !!---------------------------------------------------------------------- 
     22   USE oce_trc        !  
     23   USE trc            ! 
     24   USE daymod         ! 
     25   USE iom            ! 
    2526    
    2627   IMPLICIT NONE 
    2728   PRIVATE 
    2829 
    29    PUBLIC   trc_rst_opn       ! called by ??? 
    30    PUBLIC   trc_rst_read      ! called by ??? 
    31    PUBLIC   trc_rst_wri       ! called by ??? 
    32    PUBLIC   trc_rst_cal 
    33  
    34    !!---------------------------------------------------------------------- 
    35    !! NEMO/TOP 3.7 , NEMO Consortium (2018) 
     30   PUBLIC   trc_rst_opn    ! called by trcstp 
     31   PUBLIC   trc_rst_read   ! called by trcini 
     32   PUBLIC   trc_rst_wri    ! called by trcstp 
     33   PUBLIC   trc_rst_cal    ! called by trcstp & trcini (and OFF/nemogcm) 
     34 
     35   !!---------------------------------------------------------------------- 
     36   !! NEMO/TOP 4.0 , NEMO Consortium (2018) 
    3637   !! $Id$ 
    3738   !! Software governed by the CeCILL licence (./LICENSE) 
    … …  
    9293      ! 
    9394   END SUBROUTINE trc_rst_opn 
     95 
    9496 
    9597   SUBROUTINE trc_rst_read 
    … …  
    269271            ENDIF 
    270272            ! 
    271             IF( ln_rsttr )  THEN   ;    neuler = 1 
    272             ELSE                   ;    neuler = 0 
     273            IF( ln_rsttr )  THEN   ;    l_1st_euler = .FALSE.     ! OFF restart: no Euler 1st time-step 
     274            ELSE                   ;    l_1st_euler = .TRUE.      ! OFF cold start: Euler 1st time-step is used 
    273275            ENDIF 
    274276            ! 
    … …  
    346348#endif 
    347349 
    348    !!---------------------------------------------------------------------- 
    349    !! NEMO/TOP 3.3 , NEMO Consortium (2018) 
    350    !! $Id$ 
    351    !! Software governed by the CeCILL licence (./LICENSE) 
    352350   !!====================================================================== 
    353351END MODULE trcrst 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/trcstp.F90

    r9598 r9863  
    6161      IF( ln_timing )   CALL timing_start('trc_stp') 
    6262      ! 
    63       IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN     ! at nittrc000 
     63      IF( l_1st_euler .OR. ln_top_euler ) THEN     ! at nittrc000 
    6464         r2dttrc =  rdttrc           ! = rdttrc (use or restarting with Euler time stepping) 
    6565      ELSEIF( kt <= nittrc000 + nn_dttrc ) THEN          ! at nittrc000 or nittrc000+1 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/TOP/trcsub.F90

    r9598 r9863  
    466466      ! 
    467467      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    468       REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
     468      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z1_2rau0   ! local scalars 
    469469      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 
    470470      !!--------------------------------------------------------------------- 
    … …  
    486486      CALL div_hor( kt )                              ! Horizontal divergence & Relative vorticity 
    487487      ! 
    488       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    489       IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt 
    490  
    491488      !                                           !------------------------------! 
    492489      !                                           !   After Sea Surface Height   ! 
    … …  
    499496      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 
    500497      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 
    501       z1_rau0 = 0.5 / rau0 
    502       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     498      z1_2rau0 = 0.5 * r1_rau0 
     499      ssha(:,:) = (  sshb(:,:) - r2dt * ( z1_2rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    503500 
    504501      IF( .NOT.ln_dynspg_ts ) THEN 
    … …  
    517514      !                                           !     Now Vertical Velocity    ! 
    518515      !                                           !------------------------------! 
    519       z1_2dt = 1.e0 / z2dt 
    520516      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    521517         ! - ML - need 3 lines here because replacement of e3t by its expression yields too long lines otherwise 
    522518         wn(:,:,jk) = wn(:,:,jk+1) -   e3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    523519            &                      - ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )    & 
    524             &                         * tmask(:,:,jk) * z1_2dt 
     520            &                         * tmask(:,:,jk) * r1_2dt 
    525521         IF( ln_bdy ) wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
    526522      END DO 
Note: See TracChangeset for help on using the changeset viewer.