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

Changeset 6079


Ignore:
Timestamp:
2015-12-17T11:08:49+01:00 (8 years ago)
Author:
jamesharle
Message:

merge to trunk@5936

Location:
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO
Files:
5 deleted
56 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_2

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r5901 r6079  
    2222   USE oce 
    2323   USE dom_oce       
    24    USE sol_oce 
    2524   USE agrif_oce 
    2625   USE phycst 
     
    2928   USE lib_mpp 
    3029   USE wrk_nemo 
    31    USE dynspg_oce 
    3230   USE zdf_oce 
    3331  
     
    3836 
    3937   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts 
    40    PUBLIC   interpun, interpvn, interpun2d, interpvn2d  
     38   PUBLIC   interpun, interpvn 
    4139   PUBLIC   interptsn,  interpsshn 
    4240   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b 
     
    8078      !! 
    8179      INTEGER :: ji,jj,jk, j1,j2, i1,i2 
    82       REAL(wp) :: timeref 
    83       REAL(wp) :: z2dt, znugdt 
    84       REAL(wp) :: zrhox, zrhoy 
    85       REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1 
     80      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb 
    8681      !!----------------------------------------------------------------------   
    8782 
    8883      IF( Agrif_Root() )   RETURN 
    8984 
    90       CALL wrk_alloc( jpi, jpj, spgv1, spgu1 ) 
     85      CALL wrk_alloc( jpi, jpj, zub, zvb ) 
    9186 
    9287      Agrif_SpecialValue=0. 
     
    9691      CALL Agrif_Bc_variable(vn_interp_id,procname=interpvn) 
    9792 
    98 #if defined key_dynspg_flt 
    99       CALL Agrif_Bc_variable(e1u_id,calledweight=1., procname=interpun2d) 
    100       CALL Agrif_Bc_variable(e2v_id,calledweight=1., procname=interpvn2d) 
    101 #endif 
    102  
    10393      Agrif_UseSpecialValue = .FALSE. 
    104  
    105       zrhox = Agrif_Rhox() 
    106       zrhoy = Agrif_Rhoy() 
    107  
    108       timeref = 1. 
    109       ! time step: leap-frog 
    110       z2dt = 2. * rdt 
    111       ! time step: Euler if restart from rest 
    112       IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 
    113       ! coefficients 
    114       znugdt =  grav * z2dt     
    115  
     94  
    11695      ! prevent smoothing in ghost cells 
    11796      i1=1 
     
    126105 
    127106      IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    128 #if defined key_dynspg_flt 
    129          DO jk=1,jpkm1 
     107 
     108         ! Smoothing 
     109         ! --------- 
     110         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 
     111            ua_b(2,:)=0._wp 
     112            DO jk=1,jpkm1 
     113               DO jj=1,jpj 
     114                  ua_b(2,jj) = ua_b(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk) 
     115               END DO 
     116            END DO 
     117            DO jj=1,jpj 
     118               ua_b(2,jj) = ua_b(2,jj) * hur_a(2,jj)             
     119            END DO 
     120         ENDIF 
     121 
     122         DO jk=1,jpkm1                 ! Smooth 
    130123            DO jj=j1,j2 
    131                ua(2,jj,jk) = (ua(2,jj,jk) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,jk) 
    132             END DO 
    133          END DO 
    134  
    135          spgu(2,:)=0. 
    136  
     124               ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk)) 
     125               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk) 
     126            END DO 
     127         END DO 
     128 
     129         zub(2,:)=0._wp                ! Correct transport 
    137130         DO jk=1,jpkm1 
    138131            DO jj=1,jpj 
    139                spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk) 
    140             END DO 
    141          END DO 
    142  
     132               zub(2,jj) = zub(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk) 
     133            END DO 
     134         END DO 
    143135         DO jj=1,jpj 
    144             IF (umask(2,jj,1).NE.0.) THEN 
    145                spgu(2,jj)=spgu(2,jj)/hu(2,jj) 
    146             ENDIF 
    147          END DO 
    148 #else 
    149          spgu(2,:) = ua_b(2,:) 
    150 #endif 
    151  
    152          DO jk=1,jpkm1 
     136            zub(2,jj) = zub(2,jj) * hur_a(2,jj) 
     137         END DO 
     138 
     139         DO jk=1,jpkm1 
     140            DO jj=1,jpj 
     141               ua(2,jj,jk) = (ua(2,jj,jk)+ua_b(2,jj)-zub(2,jj))*umask(2,jj,jk) 
     142            END DO 
     143         END DO 
     144 
     145         ! Set tangential velocities to time splitting estimate 
     146         !----------------------------------------------------- 
     147         IF ( ln_dynspg_ts) THEN 
     148            zvb(2,:)=0._wp 
     149            DO jk=1,jpkm1 
     150               DO jj=1,jpj 
     151                  zvb(2,jj) = zvb(2,jj) + fse3v_a(2,jj,jk) * va(2,jj,jk) 
     152               END DO 
     153            END DO 
     154            DO jj=1,jpj 
     155               zvb(2,jj) = zvb(2,jj) * hvr_a(2,jj) 
     156            END DO 
     157            DO jk=1,jpkm1 
     158               DO jj=1,jpj 
     159                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj))*vmask(2,jj,jk) 
     160               END DO 
     161            END DO 
     162         ENDIF 
     163 
     164         ! Mask domain edges: 
     165         !------------------- 
     166         DO jk=1,jpkm1 
     167            DO jj=1,jpj 
     168               ua(1,jj,jk) = 0._wp 
     169               va(1,jj,jk) = 0._wp 
     170            END DO 
     171         END DO          
     172 
     173      ENDIF 
     174 
     175      IF((nbondi == 1).OR.(nbondi == 2)) THEN 
     176 
     177         ! Smoothing 
     178         ! --------- 
     179         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 
     180            ua_b(nlci-2,:)=0._wp 
     181            DO jk=1,jpkm1 
     182               DO jj=1,jpj 
     183                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk) 
     184               END DO 
     185            END DO 
     186            DO jj=1,jpj 
     187               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * hur_a(nlci-2,jj)             
     188            END DO 
     189         ENDIF 
     190 
     191         DO jk=1,jpkm1                 ! Smooth 
    153192            DO jj=j1,j2 
    154                ua(2,jj,jk) = 0.25*(ua(1,jj,jk)+2.*ua(2,jj,jk)+ua(3,jj,jk)) 
    155                ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk) 
    156             END DO 
    157          END DO 
    158  
    159          spgu1(2,:)=0. 
    160  
     193               ua(nlci-2,jj,jk) = 0.25_wp*(ua(nlci-3,jj,jk)+2._wp*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk)) 
     194               ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk) 
     195            END DO 
     196         END DO 
     197 
     198         zub(nlci-2,:)=0._wp           ! Correct transport 
    161199         DO jk=1,jpkm1 
    162200            DO jj=1,jpj 
    163                spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk) 
    164             END DO 
    165          END DO 
    166  
     201               zub(nlci-2,jj) = zub(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk) 
     202            END DO 
     203         END DO 
    167204         DO jj=1,jpj 
    168             IF (umask(2,jj,1).NE.0.) THEN 
    169                spgu1(2,jj)=spgu1(2,jj)/hu(2,jj) 
    170             ENDIF 
    171          END DO 
    172  
    173          DO jk=1,jpkm1 
    174             DO jj=j1,j2 
    175                ua(2,jj,jk) = (ua(2,jj,jk)+spgu(2,jj)-spgu1(2,jj))*umask(2,jj,jk) 
    176             END DO 
    177          END DO 
    178  
    179 #if defined key_dynspg_ts 
     205            zub(nlci-2,jj) = zub(nlci-2,jj) * hur_a(nlci-2,jj) 
     206         END DO 
     207 
     208         DO jk=1,jpkm1 
     209            DO jj=1,jpj 
     210               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+ua_b(nlci-2,jj)-zub(nlci-2,jj))*umask(nlci-2,jj,jk) 
     211            END DO 
     212         END DO 
     213 
    180214         ! Set tangential velocities to time splitting estimate 
    181          spgv1(2,:)=0. 
    182          DO jk=1,jpkm1 
     215         !----------------------------------------------------- 
     216         IF ( ln_dynspg_ts) THEN 
     217            zvb(nlci-1,:)=0._wp 
     218            DO jk=1,jpkm1 
     219               DO jj=1,jpj 
     220                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + fse3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk) 
     221               END DO 
     222            END DO 
    183223            DO jj=1,jpj 
    184                spgv1(2,jj)=spgv1(2,jj)+fse3v_a(2,jj,jk)*va(2,jj,jk) 
    185             END DO 
    186          END DO 
    187          DO jj=1,jpj 
    188             spgv1(2,jj)=spgv1(2,jj)*hvr_a(2,jj) 
    189          END DO 
     224               zvb(nlci-1,jj) = zvb(nlci-1,jj) * hvr_a(nlci-1,jj) 
     225            END DO 
     226            DO jk=1,jpkm1 
     227               DO jj=1,jpj 
     228                  va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-zvb(nlci-1,jj))*vmask(nlci-1,jj,jk) 
     229               END DO 
     230            END DO 
     231         ENDIF 
     232 
     233         ! Mask domain edges: 
     234         !------------------- 
    190235         DO jk=1,jpkm1 
    191236            DO jj=1,jpj 
    192                va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-spgv1(2,jj))*vmask(2,jj,jk) 
    193             END DO 
    194          END DO 
    195 #endif 
    196  
    197       ENDIF 
    198  
    199       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    200 #if defined key_dynspg_flt 
    201          DO jk=1,jpkm1 
    202             DO jj=j1,j2 
    203                ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)- z2dt * znugdt * laplacu(nlci-2,jj))*umask(nlci-2,jj,jk) 
    204             END DO 
    205          END DO 
    206          spgu(nlci-2,:)=0. 
    207          DO jk=1,jpkm1 
    208             DO jj=1,jpj 
    209                spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk) 
    210             ENDDO 
    211          ENDDO 
    212          DO jj=1,jpj 
    213             IF (umask(nlci-2,jj,1).NE.0.) THEN 
    214                spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj) 
    215             ENDIF 
    216          END DO 
    217 #else 
    218          spgu(nlci-2,:) = ua_b(nlci-2,:) 
    219 #endif 
    220          DO jk=1,jpkm1 
    221             DO jj=j1,j2 
    222                ua(nlci-2,jj,jk) = 0.25*(ua(nlci-3,jj,jk)+2.*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk)) 
    223  
    224                ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk) 
    225  
    226             END DO 
    227          END DO 
    228          spgu1(nlci-2,:)=0. 
    229          DO jk=1,jpkm1 
    230             DO jj=1,jpj 
    231                spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk) 
    232             END DO 
    233          END DO 
    234          DO jj=1,jpj 
    235             IF (umask(nlci-2,jj,1).NE.0.) THEN 
    236                spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj) 
    237             ENDIF 
    238          END DO 
    239          DO jk=1,jpkm1 
    240             DO jj=j1,j2 
    241                ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+spgu(nlci-2,jj)-spgu1(nlci-2,jj))*umask(nlci-2,jj,jk) 
    242             END DO 
    243          END DO 
    244  
    245 #if defined key_dynspg_ts 
     237               ua(nlci-1,jj,jk) = 0._wp 
     238               va(nlci  ,jj,jk) = 0._wp 
     239            END DO 
     240         END DO  
     241 
     242      ENDIF 
     243 
     244      IF((nbondj == -1).OR.(nbondj == 2)) THEN 
     245 
     246         ! Smoothing 
     247         ! --------- 
     248         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 
     249            va_b(:,2)=0._wp 
     250            DO jk=1,jpkm1 
     251               DO ji=1,jpi 
     252                  va_b(ji,2) = va_b(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk) 
     253               END DO 
     254            END DO 
     255            DO ji=1,jpi 
     256               va_b(ji,2) = va_b(ji,2) * hvr_a(ji,2)             
     257            END DO 
     258         ENDIF 
     259 
     260         DO jk=1,jpkm1                 ! Smooth 
     261            DO ji=i1,i2 
     262               va(ji,2,jk)=0.25_wp*(va(ji,1,jk)+2._wp*va(ji,2,jk)+va(ji,3,jk)) 
     263               va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk) 
     264            END DO 
     265         END DO 
     266 
     267         zvb(:,2)=0._wp                ! Correct transport 
     268         DO jk=1,jpkm1 
     269            DO ji=1,jpi 
     270               zvb(ji,2) = zvb(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk) 
     271            END DO 
     272         END DO 
     273         DO ji=1,jpi 
     274            zvb(ji,2) = zvb(ji,2) * hvr_a(ji,2) 
     275         END DO 
     276         DO jk=1,jpkm1 
     277            DO ji=1,jpi 
     278               va(ji,2,jk) = (va(ji,2,jk)+va_b(ji,2)-zvb(ji,2))*vmask(ji,2,jk) 
     279            END DO 
     280         END DO 
     281 
    246282         ! Set tangential velocities to time splitting estimate 
    247          spgv1(nlci-1,:)=0._wp 
    248          DO jk=1,jpkm1 
    249             DO jj=1,jpj 
    250                spgv1(nlci-1,jj)=spgv1(nlci-1,jj)+fse3v_a(nlci-1,jj,jk)*va(nlci-1,jj,jk)*vmask(nlci-1,jj,jk) 
    251             END DO 
    252          END DO 
    253  
    254          DO jj=1,jpj 
    255             spgv1(nlci-1,jj)=spgv1(nlci-1,jj)*hvr_a(nlci-1,jj) 
    256          END DO 
    257  
    258          DO jk=1,jpkm1 
    259             DO jj=1,jpj 
    260                va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-spgv1(nlci-1,jj))*vmask(nlci-1,jj,jk) 
    261             END DO 
    262          END DO 
    263 #endif 
    264  
    265       ENDIF 
    266  
    267       IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    268  
    269 #if defined key_dynspg_flt 
    270          DO jk=1,jpkm1 
     283         !----------------------------------------------------- 
     284         IF ( ln_dynspg_ts ) THEN 
     285            zub(:,2)=0._wp 
     286            DO jk=1,jpkm1 
     287               DO ji=1,jpi 
     288                  zub(ji,2) = zub(ji,2) + fse3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk) 
     289               END DO 
     290            END DO 
    271291            DO ji=1,jpi 
    272                va(ji,2,jk) = (va(ji,2,jk) - z2dt * znugdt * laplacv(ji,2))*vmask(ji,2,jk) 
    273             END DO 
    274          END DO 
    275  
    276          spgv(:,2)=0. 
    277  
     292               zub(ji,2) = zub(ji,2) * hur_a(ji,2) 
     293            END DO 
     294 
     295            DO jk=1,jpkm1 
     296               DO ji=1,jpi 
     297                  ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-zub(ji,2))*umask(ji,2,jk) 
     298               END DO 
     299            END DO 
     300         ENDIF 
     301 
     302         ! Mask domain edges: 
     303         !------------------- 
    278304         DO jk=1,jpkm1 
    279305            DO ji=1,jpi 
    280                spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk) 
    281             END DO 
    282          END DO 
    283  
     306               ua(ji,1,jk) = 0._wp 
     307               va(ji,1,jk) = 0._wp 
     308            END DO 
     309         END DO  
     310 
     311      ENDIF 
     312 
     313      IF((nbondj == 1).OR.(nbondj == 2)) THEN 
     314         ! Smoothing 
     315         ! --------- 
     316         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport 
     317            va_b(:,nlcj-2)=0._wp 
     318            DO jk=1,jpkm1 
     319               DO ji=1,jpi 
     320                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) 
     321               END DO 
     322            END DO 
     323            DO ji=1,jpi 
     324               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * hvr_a(ji,nlcj-2)             
     325            END DO 
     326         ENDIF 
     327 
     328         DO jk=1,jpkm1                 ! Smooth 
     329            DO ji=i1,i2 
     330               va(ji,nlcj-2,jk)=0.25_wp*(va(ji,nlcj-3,jk)+2._wp*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk)) 
     331               va(ji,nlcj-2,jk)=va(ji,nlcj-2,jk)*vmask(ji,nlcj-2,jk) 
     332            END DO 
     333         END DO 
     334 
     335         zvb(:,nlcj-2)=0._wp           ! Correct transport 
     336         DO jk=1,jpkm1 
     337            DO ji=1,jpi 
     338               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 
     339            END DO 
     340         END DO 
    284341         DO ji=1,jpi 
    285             IF (vmask(ji,2,1).NE.0.) THEN 
    286                spgv(ji,2)=spgv(ji,2)/hv(ji,2) 
    287             ENDIF 
    288          END DO 
    289 #else 
    290          spgv(:,2)=va_b(:,2) 
    291 #endif 
    292  
    293          DO jk=1,jpkm1 
    294             DO ji=i1,i2 
    295                va(ji,2,jk)=0.25*(va(ji,1,jk)+2.*va(ji,2,jk)+va(ji,3,jk)) 
    296                va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk) 
    297             END DO 
    298          END DO 
    299  
    300          spgv1(:,2)=0. 
    301  
     342            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * hvr_a(ji,nlcj-2) 
     343         END DO 
    302344         DO jk=1,jpkm1 
    303345            DO ji=1,jpi 
    304                spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk) 
    305             END DO 
    306          END DO 
    307  
    308          DO ji=1,jpi 
    309             IF (vmask(ji,2,1).NE.0.) THEN 
    310                spgv1(ji,2)=spgv1(ji,2)/hv(ji,2) 
    311             ENDIF 
    312          END DO 
    313  
    314          DO jk=1,jpkm1 
     346               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+va_b(ji,nlcj-2)-zvb(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 
     347            END DO 
     348         END DO 
     349 
     350         ! Set tangential velocities to time splitting estimate 
     351         !----------------------------------------------------- 
     352         IF ( ln_dynspg_ts ) THEN 
     353            zub(:,nlcj-1)=0._wp 
     354            DO jk=1,jpkm1 
     355               DO ji=1,jpi 
     356                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + fse3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk) 
     357               END DO 
     358            END DO 
    315359            DO ji=1,jpi 
    316                va(ji,2,jk) = (va(ji,2,jk)+spgv(ji,2)-spgv1(ji,2))*vmask(ji,2,jk) 
    317             END DO 
    318          END DO 
    319  
    320 #if defined key_dynspg_ts 
    321          ! Set tangential velocities to time splitting estimate 
    322          spgu1(:,2)=0._wp 
     360               zub(ji,nlcj-1) = zub(ji,nlcj-1) * hur_a(ji,nlcj-1) 
     361            END DO 
     362 
     363            DO jk=1,jpkm1 
     364               DO ji=1,jpi 
     365                  ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-zub(ji,nlcj-1))*umask(ji,nlcj-1,jk) 
     366               END DO 
     367            END DO 
     368         ENDIF 
     369 
     370         ! Mask domain edges: 
     371         !------------------- 
    323372         DO jk=1,jpkm1 
    324373            DO ji=1,jpi 
    325                spgu1(ji,2)=spgu1(ji,2)+fse3u_a(ji,2,jk)*ua(ji,2,jk)*umask(ji,2,jk) 
    326             END DO 
    327          END DO 
    328  
    329          DO ji=1,jpi 
    330             spgu1(ji,2)=spgu1(ji,2)*hur_a(ji,2) 
    331          END DO 
    332  
    333          DO jk=1,jpkm1 
    334             DO ji=1,jpi 
    335                ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-spgu1(ji,2))*umask(ji,2,jk) 
    336             END DO 
    337          END DO 
    338 #endif 
    339       ENDIF 
    340  
    341       IF((nbondj == 1).OR.(nbondj == 2)) THEN 
    342  
    343 #if defined key_dynspg_flt 
    344          DO jk=1,jpkm1 
    345             DO ji=1,jpi 
    346                va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)-z2dt * znugdt * laplacv(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 
    347             END DO 
    348          END DO 
    349  
    350  
    351          spgv(:,nlcj-2)=0. 
    352  
    353          DO jk=1,jpkm1 
    354             DO ji=1,jpi 
    355                spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk) 
    356             END DO 
    357          END DO 
    358  
    359          DO ji=1,jpi 
    360             IF (vmask(ji,nlcj-2,1).NE.0.) THEN 
    361                spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2) 
    362             ENDIF 
    363          END DO 
    364  
    365 #else 
    366          spgv(:,nlcj-2)=va_b(:,nlcj-2) 
    367 #endif 
    368  
    369          DO jk=1,jpkm1 
    370             DO ji=i1,i2 
    371                va(ji,nlcj-2,jk)=0.25*(va(ji,nlcj-3,jk)+2.*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk)) 
    372                va(ji,nlcj-2,jk) = va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 
    373             END DO 
    374          END DO 
    375  
    376          spgv1(:,nlcj-2)=0. 
    377  
    378          DO jk=1,jpkm1 
    379             DO ji=1,jpi 
    380                spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk) 
    381             END DO 
    382          END DO 
    383  
    384          DO ji=1,jpi 
    385             IF (vmask(ji,nlcj-2,1).NE.0.) THEN 
    386                spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2) 
    387             ENDIF 
    388          END DO 
    389  
    390          DO jk=1,jpkm1 
    391             DO ji=1,jpi 
    392                va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+spgv(ji,nlcj-2)-spgv1(ji,nlcj-2))*vmask(ji,nlcj-2,jk) 
    393             END DO 
    394          END DO 
    395  
    396 #if defined key_dynspg_ts 
    397          ! Set tangential velocities to time splitting estimate 
    398          spgu1(:,nlcj-1)=0._wp 
    399          DO jk=1,jpkm1 
    400             DO ji=1,jpi 
    401                spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)+fse3u_a(ji,nlcj-1,jk)*ua(ji,nlcj-1,jk) 
    402             END DO 
    403          END DO 
    404  
    405          DO ji=1,jpi 
    406             spgu1(ji,nlcj-1)=spgu1(ji,nlcj-1)*hur_a(ji,nlcj-1) 
    407          END DO 
    408  
    409          DO jk=1,jpkm1 
    410             DO ji=1,jpi 
    411                ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-spgu1(ji,nlcj-1))*umask(ji,nlcj-1,jk) 
    412             END DO 
    413          END DO 
    414 #endif 
    415  
    416       ENDIF 
    417       ! 
    418       CALL wrk_dealloc( jpi, jpj, spgv1, spgu1 ) 
     374               ua(ji,nlcj  ,jk) = 0._wp 
     375               va(ji,nlcj-1,jk) = 0._wp 
     376            END DO 
     377         END DO  
     378 
     379      ENDIF 
     380      ! 
     381      CALL wrk_dealloc( jpi, jpj, zub, zvb ) 
    419382      ! 
    420383   END SUBROUTINE Agrif_dyn 
     
    687650                  END DO 
    688651               END DO 
     652               tsa(nlci,j1:j2,k1:k2,jn) = 0._wp 
    689653            ENDDO 
    690654         ENDIF 
     
    706670                  END DO 
    707671               END DO 
     672               tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp 
    708673            ENDDO 
    709674         ENDIF 
     
    724689                  END DO 
    725690               END DO 
     691               tsa(1,j1:j2,k1:k2,jn) = 0._wp 
    726692            END DO 
    727693         ENDIF 
     
    742708                  END DO 
    743709               END DO 
     710               tsa(i1:i2,1,k1:k2,jn) = 0._wp 
    744711            ENDDO 
    745712         ENDIF 
     
    828795   END SUBROUTINE interpun 
    829796 
    830  
    831    SUBROUTINE interpun2d(ptab,i1,i2,j1,j2,before) 
    832       !!--------------------------------------------- 
    833       !!   *** ROUTINE interpun *** 
    834       !!---------------------------------------------     
    835       ! 
    836       INTEGER, INTENT(in) :: i1,i2,j1,j2 
    837       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 
    838       LOGICAL, INTENT(in) :: before 
    839       ! 
    840       INTEGER :: ji,jj 
    841       REAL(wp) :: ztref 
    842       REAL(wp) :: zrhoy  
    843       !!---------------------------------------------     
    844       ! 
    845       ztref = 1. 
    846  
    847       IF (before) THEN  
    848          DO jj=j1,j2 
    849             DO ji=i1,MIN(i2,nlci-1) 
    850                ptab(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj))  
    851             END DO 
    852          END DO 
    853       ELSE 
    854          zrhoy = Agrif_Rhoy() 
    855          DO jj=j1,j2 
    856             laplacu(i1:i2,jj) = ztref * (ptab(i1:i2,jj)/(zrhoy*e2u(i1:i2,jj))) !*umask(i1:i2,jj,1) 
    857          END DO 
    858       ENDIF 
    859       !  
    860    END SUBROUTINE interpun2d 
    861  
    862  
    863797   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before) 
    864798      !!--------------------------------------------- 
     
    895829      !         
    896830   END SUBROUTINE interpvn 
    897  
    898    SUBROUTINE interpvn2d(ptab,i1,i2,j1,j2,before) 
    899       !!--------------------------------------------- 
    900       !!   *** ROUTINE interpvn *** 
    901       !!---------------------------------------------     
    902       ! 
    903       INTEGER, INTENT(in) :: i1,i2,j1,j2 
    904       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 
    905       LOGICAL, INTENT(in) :: before 
    906       ! 
    907       INTEGER :: ji,jj 
    908       REAL(wp) :: zrhox  
    909       REAL(wp) :: ztref 
    910       !!---------------------------------------------     
    911       !  
    912       ztref = 1.     
    913       IF (before) THEN  
    914          !interpv entre 1 et k2 et interpv2d en jpkp1 
    915          DO jj=j1,MIN(j2,nlcj-1) 
    916             DO ji=i1,i2 
    917                ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1) 
    918             END DO 
    919          END DO 
    920       ELSE            
    921          zrhox = Agrif_Rhox() 
    922          DO ji=i1,i2 
    923             laplacv(ji,j1:j2) = ztref * (ptab(ji,j1:j2)/(zrhox*e1v(ji,j1:j2))) 
    924          END DO 
    925       ENDIF 
    926       !       
    927    END SUBROUTINE interpvn2d 
    928831 
    929832   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r5901 r6079  
    1111   USE lib_mpp 
    1212   USE wrk_nemo   
    13    USE dynspg_oce 
    1413   USE zdf_oce        ! vertical physics: ocean variables  
    1514 
     
    107106# endif 
    108107 
    109 # if defined key_dynspg_ts 
    110       IF (ln_bt_fw) THEN 
     108      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
    111109         ! Update time integrated transports 
    112110         IF (mod(nbcline,nbclineupdate) == 0) THEN 
     
    128126         ENDIF 
    129127      END IF 
    130 # endif 
    131128      ! 
    132129      nbcline = nbcline + 1 
     
    360357               ! 
    361358               ! Update barotropic velocities: 
    362 #if defined key_dynspg_ts 
    363                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    364                   zcorr = tabres(ji,jj) - un_b(ji,jj) 
    365                   ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
    366                END IF 
    367 #endif                
     359               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
     360                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     361                     zcorr = tabres(ji,jj) - un_b(ji,jj) 
     362                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
     363                  END IF 
     364               ENDIF              
    368365               un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1) 
    369366               !        
     
    425422               ! 
    426423               ! Update barotropic velocities: 
    427 #if defined key_dynspg_ts 
    428                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    429                   zcorr = tabres(ji,jj) - vn_b(ji,jj) 
    430                   vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
    431                END IF 
    432 #endif                
     424               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
     425                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     426                     zcorr = tabres(ji,jj) - vn_b(ji,jj) 
     427                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
     428                  END IF 
     429               ENDIF               
    433430               vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1) 
    434431               !        
     
    470467         END DO 
    471468      ELSE 
    472 #if ! defined key_dynspg_ts 
    473          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    474             DO jj=j1,j2 
    475                DO ji=i1,i2 
    476                   sshb(ji,jj) =   sshb(ji,jj) & 
    477                         & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
    478                END DO 
    479             END DO 
     469         IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
     470            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     471               DO jj=j1,j2 
     472                  DO ji=i1,i2 
     473                     sshb(ji,jj) =   sshb(ji,jj) & 
     474                           & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
     475                  END DO 
     476               END DO 
     477            ENDIF 
    480478         ENDIF 
    481 #endif 
     479 
    482480         DO jj=j1,j2 
    483481            DO ji=i1,i2 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90

    r5901 r6079  
    44   USE oce 
    55   USE dom_oce       
    6    USE sol_oce 
    76   USE agrif_oce 
    87   USE agrif_top_sponge 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r5901 r6079  
    104104   USE dom_oce 
    105105   USE nemogcm 
    106    USE sol_oce 
    107106   USE in_out_manager 
    108107   USE agrif_opa_update 
     
    172171   USE dom_oce 
    173172   USE nemogcm 
    174    USE sol_oce 
    175173   USE lib_mpp 
    176174   USE in_out_manager 
     
    210208   CALL Agrif_Bc_variable(vn_sponge_id,calledweight=1.,procname=interpvn_sponge) 
    211209 
    212 #if defined key_dynspg_ts 
    213210   Agrif_UseSpecialValue = .TRUE. 
    214211   CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 
    215212 
    216    Agrif_UseSpecialValue = ln_spc_dyn 
    217    CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 
    218    CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 
    219    CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 
    220    CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 
    221    ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0 
    222    ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0  
    223    ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0  
    224    ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0 
    225 #endif 
     213   IF ( ln_dynspg_ts ) THEN 
     214      Agrif_UseSpecialValue = ln_spc_dyn 
     215      CALL Agrif_Bc_variable(unb_id,calledweight=1.,procname=interpunb) 
     216      CALL Agrif_Bc_variable(vnb_id,calledweight=1.,procname=interpvnb) 
     217      CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 
     218      CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 
     219      ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0 
     220      ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0  
     221      ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0  
     222      ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0 
     223   ENDIF 
    226224 
    227225   Agrif_UseSpecialValue = .FALSE.  
     
    278276         ENDIF 
    279277      ENDIF 
     278 
     279      ! Check free surface scheme 
     280      IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 
     281         & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 
     282         WRITE(*,*) 'incompatible free surface scheme between grids' 
     283         WRITE(*,*) 'parent grid ln_dynspg_ts  :', Agrif_Parent(ln_dynspg_ts ) 
     284         WRITE(*,*) 'parent grid ln_dynspg_exp :', Agrif_Parent(ln_dynspg_exp) 
     285         WRITE(*,*) 'child grid  ln_dynspg_ts  :', ln_dynspg_ts 
     286         WRITE(*,*) 'child grid  ln_dynspg_exp :', ln_dynspg_exp 
     287         WRITE(*,*) 'those logicals should be identical'                   
     288         STOP 
     289      ENDIF 
     290 
    280291      ! check if masks and bathymetries match 
    281292      IF(ln_chk_bathy) THEN 
     
    319330   ! 
    320331   Agrif_UseSpecialValueInUpdate = .FALSE. 
    321    nbcline = 0 
     332   nbcline = 0  
    322333   lk_agrif_doupd = .FALSE. 
    323334   ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5901 r6079  
    107107      USE oce, ONLY:  zts    => tsa 
    108108      USE oce, ONLY:  zuslp  => ua   , zvslp  => va 
    109       USE oce, ONLY:  zwslpi => ua_sv , zwslpj => va_sv 
     109      USE zdf_oce, ONLY:  zwslpi => avmu , zwslpj => avmv  
    110110      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => rke 
    111111      ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OOO_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r5901 r6079  
    3434   USE zdfmxl             ! Mixed layer depth 
    3535   USE dom_oce, ONLY :   ndastp 
    36    USE sol_oce, ONLY :   gcx   ! Solver variables defined in memory 
    3736   USE in_out_manager     ! I/O manager 
    3837   USE iom                ! I/O module 
     
    114113            CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    115114#endif 
    116             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    117115            ! 
    118116            CALL iom_close( inum ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r5901 r6079  
    8686   ! 
    8787   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
    88    INTEGER,                   ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
     88   INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run) 
    8989   INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme 
    9090   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
     
    131131   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
    132132   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy) 
    133    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz      !: workspace for reading in global depth arrays (unstr.  bdy) 
     133   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz     !: workspace for reading in global depth arrays (unstr.  bdy) 
    134134   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    135135   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy) 
    136    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz     !: workspace for reading in global depth arrays (struct. bdy) 
     136   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz    !: workspace for reading in global depth arrays (struct. bdy) 
    137137!$AGRIF_DO_NOT_TREAT 
    138138   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r5620 r6079  
    2929   USE iom             ! IOM library 
    3030   USE in_out_manager  ! I/O logical units 
    31    USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 
    3231#if defined key_lim2 
    3332   USE ice_2 
     
    393392      END DO  ! ib_bdy 
    394393 
    395       ! bg jchanut tschanges 
    396394#if defined key_tide 
    397       ! Add tides if not split-explicit free surface else this is done in ts loop 
    398       IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
    399 #endif 
    400       ! end jchanut tschanges 
     395      IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                            
     396         DO ib_bdy = 1, nb_bdy    ! Tidal component added in ts loop 
     397            IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
     398               nblen => idx_bdy(ib_bdy)%nblen 
     399               nblenrim => idx_bdy(ib_bdy)%nblenrim 
     400               IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
     401               IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1)) 
     402               IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2)) 
     403               IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3)) 
     404            ENDIF 
     405         END DO 
     406      ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
     407         ! 
     408         CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
     409      ENDIF 
     410#endif 
    401411 
    402412      IF ( ln_apr_obc ) THEN 
     
    428438      !!                 
    429439      !!---------------------------------------------------------------------- 
    430       USE dynspg_oce, ONLY: lk_dynspg_ts 
    431440      !! 
    432441      INTEGER     ::  ib_bdy, jfld, jstart, jend, ierror  ! local indices 
     
    435444      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
    436445      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files 
     446      CHARACTER(len = 256)::   clname                           ! temporary file name 
    437447      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data 
    438448                                                                ! =F => baroclinic velocities in 3D boundary data 
     
    675685            ! sea ice 
    676686            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
    677  
    678687               ! Test for types of ice input (lim2 or lim3)  
    679                CALL iom_open ( bn_a_i%clname, inum ) 
    680                id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     688               ! Build file name to find dimensions  
     689               clname=TRIM(bn_a_i%clname) 
     690               IF( .NOT. bn_a_i%ln_clim ) THEN    
     691                                                  WRITE(clname, '(a,"_y",i4.4)' ) TRIM( bn_a_i%clname ), nyear    ! add year 
     692                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"m" ,i2.2)' ) TRIM( clname        ), nmonth   ! add month 
     693               ELSE 
     694                  IF( bn_a_i%cltype /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( bn_a_i%clname ), nmonth   ! add month 
     695               ENDIF 
     696               IF( bn_a_i%cltype == 'daily' .OR. bn_a_i%cltype(1:4) == 'week' ) & 
     697               &                                  WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname        ), nday     ! add day 
     698               ! 
     699               CALL iom_open  ( clname, inum ) 
     700               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
    681701               CALL iom_close ( inum ) 
    682                !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 
    683                !CALL iom_open ( bn_a_i%clname, inum ) 
    684                !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     702 
    685703                IF ( zndims == 4 ) THEN 
    686704                 ll_bdylim3 = .TRUE.   ! lim3 input 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r4792 r6079  
    2424   USE oce             ! ocean dynamics and tracers  
    2525   USE dom_oce         ! ocean space and time domain 
    26    USE dynspg_oce       
    2726   USE bdy_oce         ! ocean open boundary conditions 
    2827   USE bdydyn2d        ! open boundary conditions for barotropic solution 
     
    3534   PRIVATE 
    3635 
    37    PUBLIC   bdy_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or  
    38                         ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 
     36   PUBLIC   bdy_dyn    ! routine called in dyn_nxt 
    3937 
    4038#  include "domzgr_substitute.h90" 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r5620 r6079  
    2323   USE bdy_oce         ! ocean open boundary conditions 
    2424   USE bdylib          ! BDY library routines 
    25    USE dynspg_oce      ! for barotropic variables 
    2625   USE phycst          ! physical constants 
    2726   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r5620 r6079  
    3333!   USE tide_mod       ! Useless ?? 
    3434   USE fldread 
    35    USE dynspg_oce, ONLY: lk_dynspg_ts 
    3635 
    3736   IMPLICIT NONE 
     
    5453   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
    5554!$AGRIF_END_DO_NOT_TREAT 
    56    TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
     55   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    5756 
    5857   !!---------------------------------------------------------------------- 
     
    270269            ENDIF 
    271270            ! 
    272             IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 
    273                                      ! time splitting integration 
    274                ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
    275                ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
    276                ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
    277                dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 
    278                dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 
    279                dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 
    280             ENDIF 
     271            ! Allocate slow varying data in the case of time splitting: 
     272            ! Do it anyway because at this stage knowledge of free surface scheme is unknown 
     273            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
     274            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
     275            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
     276            dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 
     277            dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 
     278            dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 
    281279            ! 
    282280         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
     
    397395      !! 
    398396      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step 
    399       INTEGER,          DIMENSION(jpbgrd) :: ilen0  
     397      INTEGER, DIMENSION(jpbgrd) :: ilen0  
    400398      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts 
    401399      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices 
     
    416414      ! Absolute time from model initialization:    
    417415      IF( PRESENT(kit) ) THEN   
    418          z_arg = ( kt + (kit+0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt 
     416         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt 
    419417      ELSE                               
    420418         z_arg = ( kt + time_add ) * rdt 
     
    456454            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
    457455            ! 
    458             ! If time splitting, save data at first barotropic iteration 
    459             IF ( PRESENT(kit) ) THEN 
    460                IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 
    461                   IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 
    462                   IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 
    463                   IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 
    464  
    465                ELSE ! Initialize arrays from slow varying open boundary data:             
    466                   IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
    467                   IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
    468                   IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
    469                ENDIF 
     456            ! If time splitting, initialize arrays from slow varying open boundary data: 
     457            IF ( PRESENT(kit) ) THEN            
     458               IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
     459               IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
     460               IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
    470461            ENDIF 
    471462            ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r5901 r6079  
    1010   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
    1111   !!---------------------------------------------------------------------- 
    12 #if   defined key_bdy   &&   defined key_dynspg_flt 
     12#if   defined key_bdy 
    1313   !!---------------------------------------------------------------------- 
    14    !!   'key_bdy'            AND      unstructured open boundary conditions 
    15    !!   'key_dynspg_flt'                              filtered free surface 
     14   !!   'key_bdy'      unstructured open boundary conditions 
    1615   !!---------------------------------------------------------------------- 
    1716   USE oce             ! ocean dynamics and tracers  
     
    3029   PRIVATE 
    3130 
    32    PUBLIC bdy_vol        ! routine called by dynspg_flt.h90 
     31   PUBLIC bdy_vol       
    3332 
    3433   !! * Substitutions 
     
    4544      !!                      ***  ROUTINE bdyvol  *** 
    4645      !! 
    47       !! ** Purpose :   This routine is called in dynspg_flt to control  
    48       !!      the volume of the system. A correction velocity is calculated 
     46      !! ** Purpose :   This routine controls the volume of the system.  
     47      !!      A correction velocity is calculated 
    4948      !!      to correct the total transport through the unstructured OBC.  
    5049      !!      The total depth used is constant (H0) to be consistent with the  
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90

    r5901 r6079  
    1818#endif 
    1919   USE dyncor_c1d      ! Coriolis term (c1d case)         (dyn_cor_1d     ) 
    20    USE dynnxt_c1d      ! time-stepping                    (dyn_nxt routine) 
     20   USE dynnxt          ! time-stepping                    (dyn_nxt routine) 
    2121   USE dyndmp          ! U & V momentum damping           (dyn_dmp routine) 
    2222   USE restart         ! restart  
     
    139139                        CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis 
    140140                        CALL dyn_zdf    ( kstp )   ! vertical diffusion 
    141                         CALL dyn_nxt_c1d( kstp )   ! lateral velocity at next time step 
     141                        CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step 
    142142 
    143143      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r5620 r6079  
    1414   USE dom_oce         ! ocean space and time domain 
    1515   USE phycst 
    16    USE dynspg_oce 
    17    USE dynspg_ts 
    1816   USE daymod 
    1917   USE tide_mod 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r5901 r6079  
    3030   USE zdf_oce         ! ocean vertical physics 
    3131   USE ldftra          ! lateral physics: eddy diffusivity coef. 
    32    USE sol_oce         ! solver variables 
    3332   USE sbc_oce         ! Surface boundary condition: ocean fields 
    3433   USE sbc_ice         ! Surface boundary condition: ice fields 
     
    4645   USE iom 
    4746   USE ioipsl 
    48    USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities      
    4947 
    5048#if defined key_lim2 
     
    206204         CALL iom_put( "sbu", z2d )                ! bottom i-current 
    207205      ENDIF 
    208 #if defined key_dynspg_ts 
    209       CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
    210 #else 
    211       CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current 
    212 #endif 
     206 
     207      IF ( ln_dynspg_ts ) THEN 
     208         CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
     209      ELSE 
     210         CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current 
     211      ENDIF 
    213212       
    214213      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
     
    223222         CALL iom_put( "sbv", z2d )                ! bottom j-current 
    224223      ENDIF 
    225 #if defined key_dynspg_ts 
    226       CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current 
    227 #else 
    228       CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current 
    229 #endif 
     224 
     225      IF ( ln_dynspg_ts ) THEN 
     226         CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current 
     227      ELSE 
     228         CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current 
     229      ENDIF 
    230230 
    231231      CALL iom_put( "woce", wn )                   ! vertical velocity 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5901 r6079  
    4545   LOGICAL , PUBLIC ::   ln_crs          !: Apply grid coarsening to dynamical model output or online passive tracers 
    4646 
     47   !! Free surface parameters 
     48   !! ======================= 
     49   LOGICAL, PUBLIC :: ln_dynspg_exp      !: Explicit free surface flag 
     50   LOGICAL, PUBLIC :: ln_dynspg_ts       !: Split-Explicit free surface flag 
     51 
    4752   !! Time splitting parameters 
    4853   !! ========================= 
    4954   LOGICAL,  PUBLIC :: ln_bt_fw          !: Forward integration of barotropic sub-stepping 
    5055   LOGICAL,  PUBLIC :: ln_bt_av          !: Time averaging of barotropic variables 
    51    LOGICAL,  PUBLIC :: ln_bt_nn_auto     !: Set number of barotropic iterations automatically 
     56   LOGICAL,  PUBLIC :: ln_bt_auto        !: Set number of barotropic iterations automatically 
    5257   INTEGER,  PUBLIC :: nn_bt_flt         !: Filter choice 
    5358   INTEGER,  PUBLIC :: nn_baro           !: Number of barotropic iterations during one baroclinic step (rdt) 
    54    REAL(wp), PUBLIC :: rn_bt_cmax        !: Maximum allowed courant number (used if ln_bt_nn_auto=T) 
     59   REAL(wp), PUBLIC :: rn_bt_cmax        !: Maximum allowed courant number (used if ln_bt_auto=T) 
    5560 
    5661   !! Horizontal grid parameters for domhgr 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5901 r6079  
    2828   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2929   USE lib_mpp 
    30    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient 
    3130   USE wrk_nemo        ! Memory allocation 
    3231   USE timing          ! Timing 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5901 r6079  
    3535   USE dtauvd          ! data: U & V current             (dta_uvd routine) 
    3636   USE domvvl          ! varying vertical mesh 
    37    USE dynspg_oce      ! pressure gradient schemes 
    38    USE dynspg_flt      ! filtered free surface 
    39    USE sol_oce         ! ocean solver variables 
    4037   ! 
    4138   USE in_out_manager  ! I/O manager 
     
    133130      ENDIF 
    134131      ! 
    135       IF( lk_agrif ) THEN                  ! read free surface arrays in restart file 
    136          IF( ln_rstart ) THEN 
    137             IF( lk_dynspg_flt )  THEN      ! read or initialize the following fields 
    138                !                           ! gcx, gcxb for agrif_opa_init 
    139                IF( sol_oce_alloc()  > 0 )   CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed') 
    140                CALL flt_rst( nit000, 'READ' ) 
    141             ENDIF 
    142          ENDIF                             ! explicit case not coded yet with AGRIF 
    143       ENDIF 
    144       ! 
    145132      !  
    146133      ! Initialize "now" and "before" barotropic velocities: 
     
    196183      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)' 
    197184      ! 
    198 #if ! defined key_istate_fixed 
    199185      DO jk = 1, jpk 
    200186         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   & 
     
    202188         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 
    203189      END DO 
    204 #else 
    205       tsn(:,:,:,jp_tem) = ztem * tmask(:,:,:) 
    206       tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 
    207 #endif 
    208190      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 
    209191      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
    210  
    211192      ! 
    212193   END SUBROUTINE istate_t_s 
     
    451432      !!                 p=integral [ rau*g dz ] 
    452433      !!---------------------------------------------------------------------- 
    453       USE dynspg          ! surface pressure gradient             (dyn_spg routine) 
    454434      USE divhor          ! hor. divergence                       (div_hor routine) 
    455435      USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    456436      ! 
    457437      INTEGER ::   ji, jj, jk        ! dummy loop indices 
    458       INTEGER ::   indic             ! ??? 
    459438      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars 
    460439      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn 
     
    523502      vb(:,:,:) = vn(:,:,:) 
    524503       
    525       ! WARNING !!!!! 
    526       ! after initializing u and v, we need to calculate the initial streamfunction bsf. 
    527       ! Otherwise, only the trend will be computed and the model will blow up (inconsistency). 
    528       ! to do that, we call dyn_spg with a special trick: 
    529       ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the 
    530       ! right value assuming the velocities have been set up in one time step. 
    531       ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.) 
    532       !  sets up s false trend to calculate the barotropic streamfunction. 
    533  
    534       ua(:,:,:) = ub(:,:,:) / rdt 
    535       va(:,:,:) = vb(:,:,:) / rdt 
    536  
    537       ! calls dyn_spg. we assume euler time step, starting from rest. 
    538       indic = 0 
    539       CALL dyn_spg( nit000, indic )       ! surface pressure gradient 
    540       ! 
    541       ! the new velocity is ua*rdt 
    542       ! 
    543       CALL lbc_lnk( ua, 'U', -1. ) 
    544       CALL lbc_lnk( va, 'V', -1. ) 
    545  
    546       ub(:,:,:) = ua(:,:,:) * rdt 
    547       vb(:,:,:) = va(:,:,:) * rdt 
    548       ua(:,:,:) = 0.e0 
    549       va(:,:,:) = 0.e0 
    550       un(:,:,:) = ub(:,:,:) 
    551       vn(:,:,:) = vb(:,:,:) 
    552504      ! 
    553505!!gm  Check  here call to div_hor should not be necessary 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5901 r6079  
    3636   USE trd_oce         ! trends: ocean variables 
    3737   USE trddyn          ! trend manager: dynamics 
     38!jc   USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
    3839   ! 
    3940   USE in_out_manager  ! I/O manager 
     
    5859   LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme) 
    5960   LOGICAL , PUBLIC ::   ln_hpg_isf      !: s-coordinate similar to sco modify for isf 
    60    LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag 
    6161 
    6262   INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 
     
    132132      !! 
    133133      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    134          &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 
     134         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 
    135135      !!---------------------------------------------------------------------- 
    136136      ! 
     
    155155         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    156156         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    157          WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp 
    158157      ENDIF 
    159158      ! 
     
    293292      ENDIF 
    294293 
     294      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 
     295!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 
    295296 
    296297      ! Local constant initialization 
     
    491492      ! iniitialised to 0. zhpi zhpi  
    492493      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 
     494 
     495      ! Partial steps: top & bottom before horizontal gradient 
     496!jc      CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi,   &  
     497!jc               &       rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     498!jc               &      grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  ) 
    493499 
    494500!==================================================================================      
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r5901 r6079  
    1919   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    2020   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
     21   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification 
    2122   !!------------------------------------------------------------------------- 
    2223   
     
    2829   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2930   USE phycst          ! physical constants 
    30    USE dynspg_oce      ! type of surface pressure gradient 
    3131   USE dynadv          ! dynamics: vector invariant versus flux form 
    3232   USE domvvl          ! variable volume 
     
    6868      !!                  ***  ROUTINE dyn_nxt  *** 
    6969      !!                    
    70       !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary  
    71       !!             condition on the after velocity, achieved the time stepping  
     70      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary  
     71      !!             condition on the after velocity, achieve the time stepping  
    7272      !!             by applying the Asselin filter on now fields and swapping  
    7373      !!             the fields. 
    7474      !! 
    75       !! ** Method  : * After velocity is compute using a leap-frog scheme: 
    76       !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va) 
    77       !!             Note that with flux form advection and variable volume layer 
    78       !!             (lk_vvl=T), the leap-frog is applied on thickness weighted 
    79       !!             velocity. 
    80       !!             Note also that in filtered free surface (lk_dynspg_flt=T), 
    81       !!             the time stepping has already been done in dynspg module 
     75      !! ** Method  : * Ensure after velocities transport matches time splitting 
     76      !!             estimate (ln_dynspg_ts=T) 
    8277      !! 
    8378      !!              * Apply lateral boundary conditions on after velocity  
     
    10196      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    10297      INTEGER  ::   iku, ikv     ! local integers 
    103 #if ! defined key_dynspg_flt 
    104       REAL(wp) ::   z2dt         ! temporary scalar 
    105 #endif 
    106       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars 
     98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars 
    10799      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    108100      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve 
     
    112104      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt') 
    113105      ! 
    114       CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
    115       IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve ) 
     106      IF( ln_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve ) 
     107      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     108      IF( l_trddyn     )   CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 
    116109      ! 
    117110      IF( kt == nit000 ) THEN 
     
    121114      ENDIF 
    122115 
    123 #if defined key_dynspg_flt 
    124       ! 
    125       ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine 
    126       ! ------------- 
    127  
    128       ! Update after velocity on domain lateral boundaries      (only local domain required) 
    129       ! -------------------------------------------------- 
    130       CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries 
    131       CALL lbc_lnk( va, 'V', -1. )  
    132       ! 
    133 #else 
    134  
    135 # if defined key_dynspg_exp 
    136       ! Next velocity :   Leap-frog time stepping 
    137       ! ------------- 
    138       z2dt = 2. * rdt                                 ! Euler or leap-frog time step  
    139       IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt 
    140       ! 
    141       IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     116      IF ( ln_dynspg_ts ) THEN 
     117         ! Ensure below that barotropic velocities match time splitting estimate 
     118         ! Compute actual transport and replace it with ts estimate at "after" time step 
     119         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
     120         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
     121         DO jk = 2, jpkm1 
     122            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
     123            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     124         END DO 
    142125         DO jk = 1, jpkm1 
    143             ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    144             va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    145          END DO 
    146       ELSE                                            ! applied on thickness weighted velocity 
    147          DO jk = 1, jpkm1 
    148             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
    149                &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
    150                &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
    151             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
    152                &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
    153                &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    154          END DO 
    155       ENDIF 
    156 # endif 
    157  
    158 # if defined key_dynspg_ts 
    159 !!gm IF ( lk_dynspg_ts ) THEN .... 
    160       ! Ensure below that barotropic velocities match time splitting estimate 
    161       ! Compute actual transport and replace it with ts estimate at "after" time step 
    162       zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 
    163       zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 
    164       DO jk = 2, jpkm1 
    165          zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    166          zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
    167       END DO 
    168       DO jk = 1, jpkm1 
    169          ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
    170          va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
    171       END DO 
    172  
    173       IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 
    174          ! Remove advective velocity from "now velocities"  
    175          ! prior to asselin filtering      
    176          ! In the forward case, this is done below after asselin filtering    
    177          ! so that asselin contribution is removed at the same time  
    178          DO jk = 1, jpkm1 
    179             un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
    180             vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
    181          END DO   
    182       ENDIF 
    183 !!gm ENDIF 
    184 # endif 
     126            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 
     127            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 
     128         END DO 
     129 
     130         IF ( .NOT.ln_bt_fw ) THEN 
     131            ! Remove advective velocity from "now velocities"  
     132            ! prior to asselin filtering      
     133            ! In the forward case, this is done below after asselin filtering    
     134            ! so that asselin contribution is removed at the same time  
     135            DO jk = 1, jpkm1 
     136               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 
     137               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 
     138            END DO   
     139         ENDIF 
     140      ENDIF 
    185141 
    186142      ! Update after velocity on domain lateral boundaries 
    187143      ! --------------------------------------------------       
    188       CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries 
    189       CALL lbc_lnk( va, 'V', -1. )  
    190       ! 
    191 # if defined key_bdy 
    192       !                                !* BDY open boundaries 
    193       IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt ) 
    194       IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 
    195  
    196 !!$   Do we need a call to bdy_vol here?? 
    197       ! 
    198 # endif 
    199       ! 
    200144# if defined key_agrif 
    201145      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries 
    202146# endif 
    203 #endif 
    204  
     147      ! 
     148      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries 
     149      CALL lbc_lnk( va, 'V', -1. )  
     150      ! 
     151# if defined key_bdy 
     152      !                                !* BDY open boundaries 
     153      IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 
     154      IF( lk_bdy .AND. ln_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 
     155 
     156!!$   Do we need a call to bdy_vol here?? 
     157      ! 
     158# endif 
     159      ! 
    205160      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics 
    206161         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step  
     
    259214            ! (used as a now filtered scale factor until the swap) 
    260215            ! ---------------------------------------------------- 
    261             IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     216            IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 
    262217               ! No asselin filtering on thicknesses if forward time splitting 
    263                   fse3t_b(:,:,:) = fse3t_n(:,:,:) 
     218               fse3t_b(:,:,:) = fse3t_n(:,:,:) 
    264219            ELSE 
    265220               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 
     
    336291         ENDIF 
    337292         ! 
    338          IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
     293         IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 
    339294            ! Revert "before" velocities to time split estimate 
    340295            ! Doing it here also means that asselin filter contribution is removed   
     
    400355      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   & 
    401356         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    402       !  
    403       CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva ) 
    404       IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve ) 
     357      ! 
     358      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve ) 
     359      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     360      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 
    405361      ! 
    406362      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5901 r6079  
    66   !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code 
    77   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option 
    8    !!---------------------------------------------------------------------- 
    9  
    10    !!---------------------------------------------------------------------- 
    11    !!   dyn_spg     : update the dynamics trend with the lateral diffusion 
    12    !!   dyn_spg_ctl : initialization, namelist read, and parameters control 
     8   !!            3.7  ! 2015-11  (J. Chanut) Suppression of filtered free surface  
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   dyn_spg     : update the dynamics trend with surface pressure gradient  
     13   !!   dyn_spg_init: initialization, namelist read, and parameters control 
    1314   !!---------------------------------------------------------------------- 
    1415   USE oce            ! ocean dynamics and tracers variables 
     
    1819   USE sbc_oce        ! surface boundary condition: ocean 
    1920   USE sbcapr         ! surface boundary condition: atmospheric pressure 
    20    USE dynspg_oce     ! surface pressure gradient variables 
    2121   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine) 
    2222   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine) 
    23    USE dynspg_flt     ! surface pressure gradient     (dyn_spg_flt routine) 
    24    USE dynadv         ! dynamics: vector invariant versus flux form 
    25    USE dynhpg, ONLY: ln_dynhpg_imp 
    2623   USE sbctide 
    2724   USE updtide 
     
    3229   USE in_out_manager ! I/O manager 
    3330   USE lib_mpp        ! MPP library 
    34    USE solver         ! solver initialization 
    3531   USE wrk_nemo       ! Memory Allocation 
    3632   USE timing         ! Timing 
     
    4339   PUBLIC   dyn_spg_init   ! routine called by opa module 
    4440 
    45    INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...  
     41   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from ln_dynspg_...  
    4642 
    4743   !! * Substitutions 
     
    5551CONTAINS 
    5652 
    57    SUBROUTINE dyn_spg( kt, kindic ) 
     53   SUBROUTINE dyn_spg( kt ) 
    5854      !!---------------------------------------------------------------------- 
    5955      !!                  ***  ROUTINE dyn_spg  *** 
     
    6662      !!gm            the after velocity, in the 2 other (ua,va) are still the trends 
    6763      !! 
    68       !! ** Method  :   Three schemes: 
     64      !! ** Method  :   Two schemes: 
    6965      !!              - explicit computation      : the spg is evaluated at now 
    70       !!              - filtered computation      : the Roulet & madec (2000) technique is used 
    7166      !!              - split-explicit computation: a time splitting technique is used 
    7267      !! 
     
    8075      ! 
    8176      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    82       INTEGER, INTENT(  out) ::   kindic   ! solver flag 
    8377      ! 
    8478      INTEGER  ::   ji, jj, jk                             ! dummy loop indices 
     
    10397 
    10498      IF(      ln_apr_dyn                                                &   ! atmos. pressure 
    105          .OR.  ( .NOT.lk_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
     99         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) )   &   ! tide potential (no time slitting) 
    106100         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice 
    107101         ! 
     
    113107         END DO          
    114108         ! 
    115          IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==! 
     109         IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN                    !==  Atmospheric pressure gradient (added later in time-split case) ==! 
    116110            zg_2 = grav * 0.5 
    117111            DO jj = 2, jpjm1                          ! gradient of Patm using inverse barometer ssh 
     
    126120         ! 
    127121         !                                    !==  tide potential forcing term  ==! 
    128          IF( .NOT.lk_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
     122         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case 
    129123            ! 
    130124            CALL upd_tide( kt )                      ! update tide potential 
     
    171165      CASE (  0 )   ;   CALL dyn_spg_exp( kt )              ! explicit 
    172166      CASE (  1 )   ;   CALL dyn_spg_ts ( kt )              ! time-splitting 
    173       CASE (  2 )   ;   CALL dyn_spg_flt( kt, kindic )      ! filtered 
    174167      !                                                     
    175168      END SELECT 
    176169      !                     
    177170      IF( l_trddyn )   THEN                      ! save the surface pressure gradient trends for further diagnostics 
    178          SELECT CASE ( nspg ) 
    179          CASE ( 0, 1 ) 
    180             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    181             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    182          CASE( 2 ) 
    183             z2dt = 2. * rdt 
    184             IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    185             ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 
    186             ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 
    187          END SELECT 
     171         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     172         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    188173         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 
    189174         ! 
     
    203188      !!                  ***  ROUTINE dyn_spg_init  *** 
    204189      !!                 
    205       !! ** Purpose :   Control the consistency between cpp options for  
     190      !! ** Purpose :   Control the consistency between namelist options for  
    206191      !!              surface pressure gradient schemes 
    207192      !!---------------------------------------------------------------------- 
    208       INTEGER ::   ioptio 
     193      INTEGER ::   ioptio, ios 
     194      ! 
     195      NAMELIST/namdyn_spg/ ln_dynspg_exp, ln_dynspg_ts, & 
     196      &                    ln_bt_fw, ln_bt_av, ln_bt_auto, & 
     197      &                    nn_baro, rn_bt_cmax, nn_bt_flt 
    209198      !!---------------------------------------------------------------------- 
    210199      ! 
    211200      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_init') 
    212201      ! 
    213       IF(lwp) THEN             ! Control print 
     202      REWIND( numnam_ref )              ! Namelist namdyn_spg in reference namelist : Free surface 
     203      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 
     204901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp ) 
     205 
     206      REWIND( numnam_cfg )              ! Namelist namdyn_spg in configuration namelist : Free surface 
     207      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 
     208902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 
     209      IF(lwm) WRITE ( numond, namdyn_spg ) 
     210      ! 
     211      IF(lwp) THEN             ! Namelist print 
    214212         WRITE(numout,*) 
    215213         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 
    216214         WRITE(numout,*) '~~~~~~~~~~~' 
    217          WRITE(numout,*) '     Explicit free surface                  lk_dynspg_exp = ', lk_dynspg_exp 
    218          WRITE(numout,*) '     Free surface with time splitting       lk_dynspg_ts  = ', lk_dynspg_ts 
    219          WRITE(numout,*) '     Filtered free surface cst volume       lk_dynspg_flt = ', lk_dynspg_flt 
    220       ENDIF 
    221  
    222       IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 
    223       ! (do it now, to set nn_baro, used to allocate some arrays later on) 
    224       !                        ! allocate dyn_spg arrays 
    225       IF( lk_dynspg_ts ) THEN 
    226          IF( dynspg_oce_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays') 
     215         WRITE(numout,*) '     Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp 
     216         WRITE(numout,*) '     Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts 
     217      ENDIF 
     218 
     219      IF( ln_dynspg_ts ) THEN 
     220         CALL dyn_spg_ts_init( nit000  ) ! do it first, to set nn_baro, used to allocate some arrays later on 
     221         !                               ! allocate dyn_spg arrays 
    227222         IF( dyn_spg_ts_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays') 
    228223         IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 
     
    231226      !                        ! Control of surface pressure gradient scheme options 
    232227      ioptio = 0 
    233       IF(lk_dynspg_exp)   ioptio = ioptio + 1 
    234       IF(lk_dynspg_ts )   ioptio = ioptio + 1 
    235       IF(lk_dynspg_flt)   ioptio = ioptio + 1 
     228      IF(ln_dynspg_exp)   ioptio = ioptio + 1 
     229      IF(ln_dynspg_ts )   ioptio = ioptio + 1 
    236230      ! 
    237231      IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    238            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 
    239       IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   & 
    240            &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 
    241       ! 
    242       IF( lk_dynspg_exp)   nspg =  0 
    243       IF( lk_dynspg_ts )   nspg =  1 
    244       IF( lk_dynspg_flt)   nspg =  2 
     232           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 
     233      IF( ln_dynspg_ts .AND. ln_isfcav )   & 
     234           &   CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 
     235      ! 
     236      IF( ln_dynspg_exp)   nspg =  0 
     237      IF( ln_dynspg_ts )   nspg =  1 
    245238      ! 
    246239      IF(lwp) THEN 
     
    248241         IF( nspg ==  0 )   WRITE(numout,*) '     explicit free surface' 
    249242         IF( nspg ==  1 )   WRITE(numout,*) '     free surface with time splitting scheme' 
    250          IF( nspg ==  2 )   WRITE(numout,*) '     filtered free surface' 
    251       ENDIF 
    252  
    253 #if defined key_dynspg_flt 
    254       CALL solver_init( nit000 )   ! Elliptic solver initialisation 
    255 #endif 
    256       !               ! Control of hydrostatic pressure choice 
    257       IF( lk_dynspg_ts .AND. ln_dynhpg_imp )   CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 
     243      ENDIF 
    258244      ! 
    259245      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init') 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r5901 r6079  
    77   !!            3.2  !  2009-06  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 
    88   !!---------------------------------------------------------------------- 
    9 #if defined key_dynspg_exp 
    109   !!---------------------------------------------------------------------- 
    11    !!   'key_dynspg_exp'                              explicit free surface 
     10   !!                     explicit free surface 
    1211   !!---------------------------------------------------------------------- 
    1312   !!   dyn_spg_exp  : update the momentum trend with the surface  
     
    3130   PRIVATE 
    3231 
    33    PUBLIC   dyn_spg_exp   ! routine called by step.F90 
     32   PUBLIC   dyn_spg_exp   ! routine called by dyn_spg  
    3433 
    3534   !! * Substitutions 
     
    101100   END SUBROUTINE dyn_spg_exp 
    102101 
    103 #else 
    104    !!---------------------------------------------------------------------- 
    105    !!   Default case :   Empty module   No standart explicit free surface  
    106    !!---------------------------------------------------------------------- 
    107 CONTAINS 
    108    SUBROUTINE dyn_spg_exp( kt )       ! Empty routine 
    109       WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 
    110    END SUBROUTINE dyn_spg_exp 
    111 #endif 
    112  
    113102   !!====================================================================== 
    114103END MODULE dynspg_exp 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5901 r6079  
    1111   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
    1212   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
     13   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
    1314   !!--------------------------------------------------------------------- 
    14 #if defined key_dynspg_ts 
    1515   !!---------------------------------------------------------------------- 
    16    !!   'key_dynspg_ts'         split explicit free surface 
     16   !!                       split explicit free surface 
    1717   !!---------------------------------------------------------------------- 
    1818   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
     
    2323   USE sbc_oce         ! surface boundary condition: ocean 
    2424   USE sbcisf          ! ice shelf variable (fwfisf) 
    25    USE dynspg_oce      ! surface pressure gradient variables 
    2625   USE phycst          ! physical constants 
    2726   USE dynvor          ! vorticity term 
    2827   USE bdy_par         ! for lk_bdy 
    29    USE bdytides        ! open boundary condition data      
     28   USE bdytides        ! open boundary condition data 
    3029   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    3130   USE sbctide         ! tides 
     
    6867   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   ! (only used with een vorticity scheme) 
    6968 
    70    ! Arrays below are saved to allow testing of the "no time averaging" option 
    71    ! If this option is not retained, these could be replaced by temporary arrays 
    72    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 
    73                                                    ubb_e, ub_e,     & 
    74                                                    vbb_e, vb_e 
    75  
    7669   !! * Substitutions 
    7770#  include "domzgr_substitute.h90" 
     
    8881      !!                  ***  routine dyn_spg_ts_alloc  *** 
    8982      !!---------------------------------------------------------------------- 
    90       INTEGER :: ierr(3) 
     83      INTEGER :: ierr(4) 
    9184      !!---------------------------------------------------------------------- 
    9285      ierr(:) = 0 
    9386 
    94       ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 
    95          &      ub_e(jpi,jpj)  , vb_e(jpi,jpj)   , & 
    96          &      ubb_e(jpi,jpj) , vbb_e(jpi,jpj)  , STAT= ierr(1) ) 
     87      ALLOCATE( ssha_e(jpi,jpj),  sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 
     88         &        ua_e(jpi,jpj),    un_e(jpi,jpj),   ub_e(jpi,jpj),   ubb_e(jpi,jpj), & 
     89         &        va_e(jpi,jpj),    vn_e(jpi,jpj),   vb_e(jpi,jpj),   vbb_e(jpi,jpj), & 
     90         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT= ierr(1) ) 
    9791 
    9892      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
     
    10195         &                          ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
    10296 
     97      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 
     98#if defined key_agrif 
     99         &      ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                              , & 
     100#endif 
     101         &      STAT= ierr(4)) 
     102 
    103103      dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
    104104 
    105105      IF( lk_mpp                )   CALL mpp_sum( dyn_spg_ts_alloc ) 
    106       IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') 
     106      IF( dyn_spg_ts_alloc /= 0 )   CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 
    107107      ! 
    108108   END FUNCTION dyn_spg_ts_alloc 
     
    146146      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
    147147      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    148       REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
    149       REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      - 
    150       REAL(wp) ::   zu_spg, zv_spg             !   -      - 
    151       REAL(wp) ::   zhura, zhvra               !   -      - 
    152       REAL(wp) ::   za0, za1, za2, za3           !   -      - 
    153       ! 
    154       REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     148      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     149      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2  !   -      - 
     150      REAL(wp) ::   zu_spg, zv_spg              !   -      - 
     151      REAL(wp) ::   zhura, zhvra          !   -      - 
     152      REAL(wp) ::   za0, za1, za2, za3    !   -      - 
     153      ! 
     154      REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 
    155155      REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 
    156       REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 
     156      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 
    157157      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 
    158158      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
     
    164164      !                                         !* Allocate temporary arrays 
    165165      CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 
    166       CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e  ) 
    167       CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 
     166      CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 
     167      CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 
    168168      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
    169169      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     
    188188      ! 
    189189                                                       ! time offset in steps for bdy data update 
    190       IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
     190      IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
    191191      ! 
    192192      IF( kt == nit000 ) THEN                !* initialisation 
     
    485485      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    486486         sshn_e(:,:) = sshn (:,:)             
    487          zun_e (:,:) = un_b (:,:)             
    488          zvn_e (:,:) = vn_b (:,:) 
     487         un_e (:,:) = un_b (:,:)             
     488         vn_e (:,:) = vn_b (:,:) 
    489489         ! 
    490490         hu_e  (:,:) = hu   (:,:)        
     
    494494      ELSE                                ! CENTRED integration: start from BEFORE fields 
    495495         sshn_e(:,:) = sshb (:,:) 
    496          zun_e (:,:) = ub_b (:,:)          
    497          zvn_e (:,:) = vb_b (:,:) 
     496         un_e (:,:) = ub_b (:,:)          
     497         vn_e (:,:) = vb_b (:,:) 
    498498         ! 
    499499         hu_e  (:,:) = hu_b (:,:)        
     
    509509      va_b  (:,:) = 0._wp 
    510510      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level 
    511       zu_sum(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    512       zv_sum(:,:) = 0._wp 
     511      un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
     512      vn_adv(:,:) = 0._wp 
    513513      !                                             ! ==================== ! 
    514514      DO jn = 1, icycle                             !  sub-time-step loop  ! 
     
    518518         ! Update only tidal forcing at open boundaries 
    519519#if defined key_tide 
    520          IF ( lk_bdy .AND. lk_tide )      CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
    521          IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset ) 
     520         IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
     521         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 
    522522#endif 
    523523         ! 
     
    534534 
    535535         ! Extrapolate barotropic velocities at step jit+0.5: 
    536          ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
    537          va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     536         ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
     537         va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
    538538 
    539539         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
     
    597597         ! Sum over sub-time-steps to compute advective velocities 
    598598         za2 = wgtbtp2(jn) 
    599          zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 
    600          zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 
     599         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 
     600         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 
    601601         ! 
    602602         ! Set next sea level: 
     
    733733         ! 
    734734         ! Add bottom stresses: 
    735          zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
    736          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     735         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
     736         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
    737737         ! 
    738738         ! Surface pressure trend: 
     
    751751            DO jj = 2, jpjm1 
    752752               DO ji = fs_2, fs_jpim1   ! vector opt. 
    753                   ua_e(ji,jj) = (                                zun_e(ji,jj)   &  
     753                  ua_e(ji,jj) = (                                 un_e(ji,jj)   &  
    754754                            &     + rdtbt * (                      zwx(ji,jj)   & 
    755755                            &                                 + zu_trd(ji,jj)   & 
     
    757757                            &   ) * umask(ji,jj,1) 
    758758 
    759                   va_e(ji,jj) = (                                zvn_e(ji,jj)   & 
     759                  va_e(ji,jj) = (                                 vn_e(ji,jj)   & 
    760760                            &     + rdtbt * (                      zwy(ji,jj)   & 
    761761                            &                                 + zv_trd(ji,jj)   & 
     
    772772                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
    773773 
    774                   ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
     774                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   &  
    775775                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
    776776                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
     
    778778                            &   ) * zhura 
    779779 
    780                   va_e(ji,jj) = (                hv_e(ji,jj)  *  zvn_e(ji,jj)   & 
     780                  va_e(ji,jj) = (                hv_e(ji,jj)  *   vn_e(ji,jj)   & 
    781781                            &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   & 
    782782                            &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   & 
     
    802802#if defined key_bdy   
    803803                                                           ! open boundaries 
    804          IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e ) 
     804         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 
    805805#endif 
    806806#if defined key_agrif                                                            
     
    810810         !                                             !  ---- 
    811811         ubb_e  (:,:) = ub_e  (:,:) 
    812          ub_e   (:,:) = zun_e (:,:) 
    813          zun_e  (:,:) = ua_e  (:,:) 
     812         ub_e   (:,:) = un_e (:,:) 
     813         un_e   (:,:) = ua_e  (:,:) 
    814814         ! 
    815815         vbb_e  (:,:) = vb_e  (:,:) 
    816          vb_e   (:,:) = zvn_e (:,:) 
    817          zvn_e  (:,:) = va_e  (:,:) 
     816         vb_e   (:,:) = vn_e (:,:) 
     817         vn_e   (:,:) = va_e  (:,:) 
    818818         ! 
    819819         sshbb_e(:,:) = sshb_e(:,:) 
     
    840840      ! ----------------------------------------------------------------------------- 
    841841      ! 
    842       ! At this stage ssha holds a time averaged value 
    843       !                                                ! Sea Surface Height at u-,v- and f-points 
    844       IF( lk_vvl ) THEN                                ! (required only in key_vvl case) 
     842      ! Set advection velocity correction: 
     843      zwx(:,:) = un_adv(:,:) 
     844      zwy(:,:) = vn_adv(:,:) 
     845      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
     846         un_adv(:,:) = zwx(:,:)*hur(:,:) 
     847         vn_adv(:,:) = zwy(:,:)*hvr(:,:) 
     848      ELSE 
     849         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 
     850         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 
     851      END IF 
     852 
     853      IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
     854         ub2_b(:,:) = zwx(:,:) 
     855         vb2_b(:,:) = zwy(:,:) 
     856      ENDIF 
     857      ! 
     858      ! Update barotropic trend: 
     859      IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
     860         DO jk=1,jpkm1 
     861            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
     862            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
     863         END DO 
     864      ELSE 
     865         ! At this stage, ssha has been corrected: compute new depths at velocity points 
    845866         DO jj = 1, jpjm1 
    846867            DO ji = 1, jpim1      ! NO Vector Opt. 
     
    854875         END DO 
    855876         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    856       ENDIF 
    857       ! 
    858       ! Set advection velocity correction: 
    859       IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
    860          un_adv(:,:) = zu_sum(:,:)*hur(:,:) 
    861          vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 
    862       ELSE 
    863          un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 
    864          vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 
    865       END IF 
    866  
    867       IF (ln_bt_fw) THEN ! Save integrated transport for next computation 
    868          ub2_b(:,:) = zu_sum(:,:) 
    869          vb2_b(:,:) = zv_sum(:,:) 
    870       ENDIF 
    871       ! 
    872       ! Update barotropic trend: 
    873       IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 
    874          DO jk=1,jpkm1 
    875             ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 
    876             va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 
    877          END DO 
    878       ELSE 
     877         ! 
    879878         DO jk=1,jpkm1 
    880879            ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 
     
    915914      ! 
    916915      CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 
    917       CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 
    918       CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 
     916      CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 
     917      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 
    919918      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 
    920919      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
     
    10641063      ! 
    10651064      INTEGER  :: ji ,jj 
    1066       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    10671065      REAL(wp) :: zxr2, zyr2, zcmax 
    10681066      REAL(wp), POINTER, DIMENSION(:,:) :: zcu 
    10691067      !! 
    1070       NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 
    1071       &                  nn_baro, rn_bt_cmax, nn_bt_flt 
    10721068      !!---------------------------------------------------------------------- 
    10731069      ! 
    1074       REWIND( numnam_ref )              ! Namelist namsplit in reference namelist : time splitting parameters 
    1075       READ  ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 
    1076 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 
    1077  
    1078       REWIND( numnam_cfg )              ! Namelist namsplit in configuration namelist : time splitting parameters 
    1079       READ  ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 
    1080 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 
    1081       IF(lwm) WRITE ( numond, namsplit ) 
    1082       ! 
    1083       !         ! Max courant number for ext. grav. waves 
     1070      ! Max courant number for ext. grav. waves 
    10841071      ! 
    10851072      CALL wrk_alloc( jpi, jpj, zcu ) 
    10861073      ! 
    1087       IF (lk_vvl) THEN  
    1088          DO jj = 1, jpj 
    1089             DO ji =1, jpi 
    1090                zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 
    1091                zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 
    1092                zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 
    1093             END DO 
    1094          END DO 
    1095       ELSE 
    1096          DO jj = 1, jpj 
    1097             DO ji =1, jpi 
    1098                zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 
    1099                zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 
    1100                zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 
    1101             END DO 
    1102          END DO 
    1103       ENDIF 
    1104  
     1074      DO jj = 1, jpj 
     1075         DO ji =1, jpi 
     1076            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 
     1077            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 
     1078            zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 
     1079         END DO 
     1080      END DO 
     1081      ! 
    11051082      zcmax = MAXVAL( zcu(:,:) ) 
    11061083      IF( lk_mpp )   CALL mpp_max( zcmax ) 
    11071084 
    11081085      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 
    1109       IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
     1086      IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
    11101087       
    11111088      rdtbt = rdt / REAL( nn_baro , wp ) 
     
    11151092      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 
    11161093      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    1117       IF( ln_bt_nn_auto ) THEN 
    1118          IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.true. Automatically set nn_baro ' 
     1094      IF( ln_bt_auto ) THEN 
     1095         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.true. Automatically set nn_baro ' 
    11191096         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax 
    11201097      ELSE 
    1121          IF(lwp) WRITE(numout,*) '     ln_ts_nn_auto=.false.: Use nn_baro in namelist ' 
     1098         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_baro in namelist ' 
    11221099      ENDIF 
    11231100 
     
    11641141   END SUBROUTINE dyn_spg_ts_init 
    11651142 
    1166 #else 
    1167    !!--------------------------------------------------------------------------- 
    1168    !!   Default case :   Empty module   No split explicit free surface 
    1169    !!--------------------------------------------------------------------------- 
    1170 CONTAINS 
    1171    INTEGER FUNCTION dyn_spg_ts_alloc()    ! Dummy function 
    1172       dyn_spg_ts_alloc = 0 
    1173    END FUNCTION dyn_spg_ts_alloc 
    1174    SUBROUTINE dyn_spg_ts( kt )            ! Empty routine 
    1175       INTEGER, INTENT(in) :: kt 
    1176       WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt 
    1177    END SUBROUTINE dyn_spg_ts 
    1178    SUBROUTINE ts_rst( kt, cdrw )          ! Empty routine 
    1179       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1180       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    1181       WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw 
    1182    END SUBROUTINE ts_rst   
    1183    SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine 
    1184       INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    1185       WRITE(*,*) 'dyn_spg_ts_init   : You should not have seen this print! error?', kt 
    1186    END SUBROUTINE dyn_spg_ts_init 
    1187 #endif 
    1188     
    11891143   !!====================================================================== 
    11901144END MODULE dynspg_ts 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r5901 r6079  
    3232   USE trd_oce        ! trends: ocean variables 
    3333   USE trddyn         ! trend manager: dynamics 
     34   USE c1d            ! 1D vertical configuration 
    3435   ! 
    3536   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    547548         END SELECT 
    548549         ! 
     550         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     551            DO jj = 1, jpjm1 
     552               DO ji = 1, fs_jpim1   ! vector opt. 
     553                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     554               END DO 
     555            END DO 
     556         ENDIF 
     557         ! 
    549558         CALL lbc_lnk( zwz, 'F', 1. ) 
    550          ! 
    551          IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
    552             DO jj = 1, jpjm1 
    553                DO ji = 1, fs_jpim1   ! vector opt. 
    554                   zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    555                END DO 
    556             END DO 
    557          ENDIF 
    558559         ! 
    559560         !                                   !==  horizontal fluxes  ==! 
     
    662663      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;    nvor_scheme = np_EEN   ;   ENDIF 
    663664      ! 
    664       IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
     665      IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 
    665666      !                       
    666667      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90

    r3625 r6079  
    88   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module 
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
     10   !!            3.7  !  2015-11  (J. Chanut) output velocities instead of trends 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    1819   USE phycst          ! physical constants 
    1920   USE zdf_oce         ! ocean vertical physics 
     21   USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 
    2022   USE sbc_oce         ! surface boundary condition: ocean 
    2123   USE lib_mpp         ! MPP library 
     
    118120         ! 
    119121      END DO                           ! End of time splitting 
     122 
     123      ! Time step momentum here to be compliant with what is done in the implicit case 
     124      ! 
     125      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     126         DO jk = 1, jpkm1 
     127            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     128            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     129         END DO 
     130      ELSE                                            ! applied on thickness weighted velocity 
     131         DO jk = 1, jpkm1 
     132            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
     133               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
     134               &         / fse3u_a(:,:,jk) * umask(:,:,jk) 
     135            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
     136               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
     137               &         / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     138         END DO 
     139      ENDIF 
    120140      ! 
    121141      CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )  
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r5901 r6079  
    2525   USE wrk_nemo        ! Memory Allocation 
    2626   USE timing          ! Timing 
    27    USE dynadv          ! dynamics: vector invariant versus flux form 
    28    USE dynspg_oce, ONLY: lk_dynspg_ts 
     27   USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 
    2928 
    3029   IMPLICIT NONE 
     
    8786      ENDIF 
    8887 
    89       ! 0. Local constant initialization 
    90       ! -------------------------------- 
    91       z1_p2dt = 1._wp / p2dt      ! inverse of the timestep 
    92  
    93       ! 1. Apply semi-implicit bottom friction 
     88      ! 1. Time step dynamics 
     89      ! --------------------- 
     90 
     91      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
     92         DO jk = 1, jpkm1 
     93            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     94            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     95         END DO 
     96      ELSE                                            ! applied on thickness weighted velocity 
     97         DO jk = 1, jpkm1 
     98            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
     99               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
     100               &                               / fse3u_a(:,:,jk) * umask(:,:,jk) 
     101            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
     102               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
     103               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk) 
     104         END DO 
     105      ENDIF 
     106 
     107      ! 2. Apply semi-implicit bottom friction 
    94108      ! -------------------------------------- 
    95109      ! Only needed for semi-implicit bottom friction setup. The explicit 
     
    97111      ! column vector of the tri-diagonal matrix equation 
    98112      ! 
    99  
    100113      IF( ln_bfrimp ) THEN 
    101114         DO jj = 2, jpjm1 
     
    119132      ENDIF 
    120133 
    121 #if defined key_dynspg_ts 
    122       IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
    123          DO jk = 1, jpkm1 
    124             ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    125             va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    126          END DO 
    127       ELSE                                            ! applied on thickness weighted velocity 
    128          DO jk = 1, jpkm1 
    129             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
    130                &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
    131                &                               / fse3u_a(:,:,jk) * umask(:,:,jk) 
    132             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
    133                &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
    134                &                               / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    135          END DO 
    136       ENDIF 
    137  
    138       IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 
     134      ! With split-explicit free surface, barotropic stress is treated explicitly 
     135      ! Update velocities at the bottom. 
     136      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
     137      !            not lead to the effective stress seen over the whole barotropic loop.  
     138      IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 
    139139         ! remove barotropic velocities: 
    140140         DO jk = 1, jpkm1 
     
    166166         END IF 
    167167      ENDIF 
    168 #endif 
    169  
    170       ! 2. Vertical diffusion on u 
     168 
     169      ! 3. Vertical diffusion on u 
    171170      ! --------------------------- 
    172171      ! Matrix and second member construction 
     
    219218      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    220219         DO ji = fs_2, fs_jpim1   ! vector opt. 
    221 #if defined key_dynspg_ts 
    222220            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
    223221            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    224222               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    225 #else 
    226             ua(ji,jj,1) = ub(ji,jj,1) & 
    227                &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    228                &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) )  
    229 #endif 
    230223         END DO 
    231224      END DO 
     
    233226         DO jj = 2, jpjm1 
    234227            DO ji = fs_2, fs_jpim1 
    235 #if defined key_dynspg_ts 
    236228               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    237 #else 
    238                zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 
    239 #endif 
    240229               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    241230            END DO 
     
    256245      END DO 
    257246 
    258 #if ! defined key_dynspg_ts 
    259       ! Normalization to obtain the general momentum trend ua 
    260       DO jk = 1, jpkm1 
    261          DO jj = 2, jpjm1    
    262             DO ji = fs_2, fs_jpim1   ! vector opt. 
    263                ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 
    264             END DO 
    265          END DO 
    266       END DO 
    267 #endif 
    268  
    269       ! 3. Vertical diffusion on v 
     247      ! 4. Vertical diffusion on v 
    270248      ! --------------------------- 
    271249      ! Matrix and second member construction 
     
    317295      ! 
    318296      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    319          DO ji = fs_2, fs_jpim1   ! vector opt. 
    320 #if defined key_dynspg_ts             
     297         DO ji = fs_2, fs_jpim1   ! vector opt.           
    321298            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
    322299            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    323300               &                                      / ( ze3va * rau0 )  
    324 #else 
    325             va(ji,jj,1) = vb(ji,jj,1) & 
    326                &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    327                &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
    328 #endif 
    329301         END DO 
    330302      END DO 
     
    332304         DO jj = 2, jpjm1 
    333305            DO ji = fs_2, fs_jpim1   ! vector opt. 
    334 #if defined key_dynspg_ts 
    335306               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    336 #else 
    337                zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 
    338 #endif 
    339307               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    340308            END DO 
     
    355323      END DO 
    356324 
    357       ! Normalization to obtain the general momentum trend va 
    358 #if ! defined key_dynspg_ts 
    359       DO jk = 1, jpkm1 
    360          DO jj = 2, jpjm1    
    361             DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt 
    363             END DO 
    364          END DO 
    365       END DO 
    366 #endif 
    367  
    368325      ! J. Chanut: Lines below are useless ? 
    369       !! restore bottom layer avmu(v)  
     326      !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 
    370327      IF( ln_bfrimp ) THEN 
    371328        DO jj = 2, jpjm1 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r5901 r6079  
    106106      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    107107 
    108 #if ! defined key_dynspg_ts 
    109       ! These lines are not necessary with time splitting since 
    110       ! boundary condition on sea level is set during ts loop 
     108      IF ( .NOT.ln_dynspg_ts ) THEN 
     109         ! These lines are not necessary with time splitting since 
     110         ! boundary condition on sea level is set during ts loop 
    111111# if defined key_agrif 
    112       CALL agrif_ssh( kt ) 
     112         CALL agrif_ssh( kt ) 
    113113# endif 
    114114# if defined key_bdy 
    115       IF( lk_bdy ) THEN 
    116          CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary 
    117          CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries 
    118       ENDIF 
     115         IF( lk_bdy ) THEN 
     116            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary 
     117            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries 
     118         ENDIF 
    119119# endif 
    120 #endif 
     120      ENDIF 
    121121 
    122122#if defined key_asminc 
     
    250250      ENDIF 
    251251 
    252 # if defined key_dynspg_ts 
    253       IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter 
    254 # else 
    255       IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter 
    256 #endif 
     252      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN  
     253                                                   !** Euler time-stepping: no filter 
    257254         sshb(:,:) = sshn(:,:)                           ! before <-- now 
    258255         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r5891 r6079  
    3232   PUBLIC   fld_map    ! routine called by tides_init 
    3333   PUBLIC   fld_read, fld_fill   ! called by sbc... modules 
     34   PUBLIC   fld_clopn 
    3435 
    3536   TYPE, PUBLIC ::   FLD_N      !: Namelist field informations 
     
    297298               ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp ) 
    298299               ztintb =  1. - ztinta 
    299 !CDIR COLLAPSE 
    300300               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2) 
    301301            ELSE   ! nothing to do... 
     
    12141214         imonth = kmonth 
    12151215         iday = kday 
     1216         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
     1217            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )   
     1218            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month 
     1219            llprevyr   = llprevmth .AND. nmonth == 1 
     1220            iyear  = nyear  - COUNT((/llprevyr /)) 
     1221            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /)) 
     1222            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday) 
     1223         ENDIF 
    12161224      ELSE                                                  ! use current day values 
    12171225         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
     
    13131321            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 
    13141322         sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
    1315          sdf(jf)%igrd     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init 
    1316          sdf(jf)%ibdy     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init 
    13171323      END DO 
    13181324 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

    r5901 r6079  
    1212   USE dom_oce         ! ocean space and time domain 
    1313   USE sbc_oce         ! surface boundary condition 
    14    USE dynspg_oce      ! surface pressure gradient variables 
    1514   USE phycst          ! physical constants 
    1615   USE fldread         ! read input fields 
     
    112111            IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data' 
    113112         ENDIF 
    114          IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts )   & 
    115             CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' ) 
     113!jc: stop below should rather be a warning  
    116114         IF( ( ln_apr_obc ) .AND. .NOT. ln_apr_dyn   )   & 
    117115            CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY requires ln_apr_dyn=T' ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90

    r5620 r6079  
    1010   USE phycst 
    1111   USE daymod 
    12    USE dynspg_oce 
    1312   USE tideini 
    1413   ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90

    r5620 r6079  
    1010   USE phycst 
    1111   USE daymod 
    12    USE dynspg_oce 
    1312   USE tide_mod 
    1413   ! 
     
    9695          &   CALL ctl_stop('rdttideramp must be positive') 
    9796       ! 
    98        IF( .NOT. lk_dynspg_ts )   CALL ctl_warn( 'sbc_tide : use of time splitting is recommended' ) 
    9997       ! 
    10098       ALLOCATE( ntide(nb_harmo) ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

    r5620 r6079  
    3131CONTAINS 
    3232 
    33    SUBROUTINE upd_tide( kt, kit, kbaro, koffset ) 
     33   SUBROUTINE upd_tide( kt, kit, time_offset ) 
    3434      !!---------------------------------------------------------------------- 
    3535      !!                 ***  ROUTINE upd_tide  *** 
     
    4242      !!----------------------------------------------------------------------       
    4343      INTEGER, INTENT(in)           ::   kt      ! ocean time-step index 
    44       INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T only) 
    45       INTEGER, INTENT(in), OPTIONAL ::   kbaro   ! number of sub-time-step           (lk_dynspg_ts=T only) 
    46       INTEGER, INTENT(in), OPTIONAL ::   koffset ! time offset in number  
    47                                                  ! of sub-time-steps                 (lk_dynspg_ts=T only) 
     44      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T) 
     45      INTEGER, INTENT(in), OPTIONAL ::   time_offset ! time offset in number  
     46                                                     ! of internal steps             (lk_dynspg_ts=F) 
     47                                                     ! of external steps             (lk_dynspg_ts=T) 
    4848      ! 
    4949      INTEGER  ::   joffset      ! local integer 
     
    5757      ! 
    5858      joffset = 0 
    59       IF( PRESENT( koffset ) )   joffset = koffset 
     59      IF( PRESENT( time_offset ) )   joffset = time_offset 
    6060      ! 
    61       IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   THEN 
    62          zt = zt + ( kit + 0.5_wp * ( joffset - 1 ) ) * rdt / REAL( kbaro, wp ) 
     61      IF( PRESENT( kit ) )   THEN 
     62         zt = zt + ( kit +  joffset - 1 ) * rdt / REAL( nn_baro, wp ) 
    6363      ELSE 
    6464         zt = zt + joffset * rdt 
     
    7474      IF( ln_tide_ramp ) THEN         ! linear increase if asked 
    7575         zt = ( kt - nit000 ) * rdt 
    76          IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp ) 
     76         IF( PRESENT( kit ) )   zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp ) 
    7777         zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  ) 
    7878         pot_astro(:,:) = zramp * pot_astro(:,:) 
     
    8686  !!---------------------------------------------------------------------- 
    8787CONTAINS 
    88   SUBROUTINE upd_tide( kt, kit, kbaro, koffset )          ! Empty routine 
     88  SUBROUTINE upd_tide( kt, kit, time_offset )  ! Empty routine 
    8989    INTEGER, INTENT(in)           ::   kt      !  integer  arg, dummy routine 
    9090    INTEGER, INTENT(in), OPTIONAL ::   kit     !  optional arg, dummy routine 
    91     INTEGER, INTENT(in), OPTIONAL ::   kbaro   !  optional arg, dummy routine 
    92     INTEGER, INTENT(in), OPTIONAL ::   koffset !  optional arg, dummy routine 
     91    INTEGER, INTENT(in), OPTIONAL ::   time_offset !  optional arg, dummy routine 
    9392    WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt 
    9493  END SUBROUTINE upd_tide 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5901 r6079  
    2626   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
    2727   USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces 
     28   USE c1d            ! 1D vertical configuration 
    2829   ! 
    2930   USE in_out_manager ! I/O manager 
     
    214215         CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 
    215216      ENDIF 
    216       IF( ioptio /= 1 )   CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 
     217      IF( (ioptio /= 1).AND. (.NOT. lk_c1d ) ) &  
     218        CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 
    217219      ! 
    218220      IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 )   &          ! Centered 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90

    r5901 r6079  
    1919   USE trd_oce        ! trends: ocean variables 
    2020   USE trdtra         ! tracers trends 
    21    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    2221   USE diaptr         ! poleward transport diagnostics 
    2322   ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90

    r5901 r6079  
    2121   USE trd_oce        ! trends: ocean variables 
    2222   USE trdtra         ! tracers trends manager 
    23    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    2423   USE sbcrnf         ! river runoffs 
    2524   USE diaptr         ! poleward transport diagnostics 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r5901 r6079  
    2020   USE trd_oce         ! trends: ocean variables 
    2121   USE trdtra          ! trends manager: tracers  
    22    USE dynspg_oce      ! surface pressure gradient variables 
    2322   USE diaptr          ! poleward transport diagnostics 
    2423   ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r5901 r6079  
    1818   USE traadv_fct      ! acces to routine interp_4th_cpt  
    1919   USE trdtra         ! trends manager: tracers  
    20    USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    2120   USE diaptr         ! poleward transport diagnostics 
    2221   ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5901 r6079  
    3131   USE zdf_oce         ! ocean vertical mixing 
    3232   USE domvvl          ! variable volume 
    33    USE dynspg_oce      ! surface     pressure gradient variables 
    34    USE dynhpg          ! hydrostatic pressure gradient  
    3533   USE trd_oce         ! trends: ocean variables 
    3634   USE trdtra          ! trends manager: tracers  
     
    5856   PUBLIC   tra_nxt_vvl   ! to be used in trcnxt 
    5957 
    60    REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg 
    6158 
    6259   !! * Substitutions 
     
    9087      !! 
    9188      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    92       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     89      !! 
    9390      !!---------------------------------------------------------------------- 
    9491      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     
    105102         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap' 
    106103         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    107          ! 
    108          rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg 
    109104      ENDIF 
    110105 
    111106      ! Update after tracer on domain lateral boundaries 
    112107      !  
     108      ! 
    113109#if defined key_agrif 
    114110      CALL Agrif_tra                     ! AGRIF zoom boundaries 
     
    181177      !!  
    182178      !! ** Method  : - Apply a Asselin time filter on now fields. 
    183       !!              - save in (ta,sa) an average over the three time levels  
    184       !!             which will be used to compute rdn and thus the semi-implicit 
    185       !!             hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    186179      !!              - swap tracer fields to prepare the next time_step. 
    187       !!                This can be summurized for tempearture as: 
    188       !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T 
    189       !!             ztm = 0                                   otherwise 
    190       !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp) 
    191       !!             tb  = tn + atfp*[ tb - 2 tn + ta ] 
    192       !!             tn  = ta   
    193       !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call) 
    194180      !! 
    195181      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    196       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     182      !!               
    197183      !!---------------------------------------------------------------------- 
    198184      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index 
     
    205191      ! 
    206192      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    207       LOGICAL  ::   ll_tra_hpg       ! local logical 
    208193      REAL(wp) ::   ztn, ztd         ! local scalars 
    209194      !!---------------------------------------------------------------------- 
     
    215200      ENDIF 
    216201      ! 
    217       IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg     
    218       ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
    219       ENDIF 
    220202      ! 
    221203      DO jn = 1, kjpt 
     
    230212                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta 
    231213                  ! 
    232                   IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average 
    233214               END DO 
    234215           END DO 
     
    248229      !!  
    249230      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
    250       !!              - save in (ta,sa) a thickness weighted average over the three  
    251       !!             time levels which will be used to compute rdn and thus the semi- 
    252       !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 
    253231      !!              - swap tracer fields to prepare the next time_step. 
    254232      !!                This can be summurized for tempearture as: 
    255       !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T 
    256       !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )    
    257       !!             ztm = 0                                                       otherwise 
    258       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 
    259       !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] ) 
    260       !!             tn  = ta  
    261       !!             ta  = zt        (NB: reset to 0 after eos_bn2 call) 
    262233      !! 
    263234      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step 
    264       !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
     235      !! 
    265236      !!---------------------------------------------------------------------- 
    266237      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
     
    276247 
    277248      !!      
    278       LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical 
     249      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical 
    279250      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
    280251      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
     
    289260      ! 
    290261      IF( cdtype == 'TRA' )  THEN    
    291          ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg 
    292262         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
    293263         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
     
    298268         END IF 
    299269      ELSE                           
    300          ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
    301270         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration 
    302271         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs 
     
    356325                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
    357326                  ! 
    358                   IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
    359                      ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
    360                      pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
    361                   ENDIF 
    362327               END DO 
    363328            END DO 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r5901 r6079  
    1818   USE zdf_oce         ! ocean vertical physics variables 
    1919   USE sbc_oce         ! surface boundary condition: ocean 
    20    USE dynspg_oce 
    2120   ! 
    2221   USE ldftra          ! lateral diffusion: eddy diffusivity 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90

    r5620 r6079  
    7171   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri = 12     !: implicit bottom friction (ln_bfrimp=.TRUE.) 
    7272   INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken  = 13     !: use for calculation of KE 
    73    INTEGER, PUBLIC, PARAMETER ::   jpdyn_spgflt  = 14  !: filter contribution to surface pressure gradient (spg_flt) 
    74    INTEGER, PUBLIC, PARAMETER ::   jpdyn_spgexp  = 15  !: explicit contribution to surface pressure gradient (spg_flt) 
    7573   ! 
    7674   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90

    r5620 r6079  
    113113      CASE( jpdyn_spg )   ;   CALL iom_put( "utrd_spg", putrd )    ! surface pressure gradient 
    114114                              CALL iom_put( "vtrd_spg", pvtrd ) 
    115       CASE( jpdyn_spgexp );   CALL iom_put( "utrd_spgexp", putrd ) ! surface pressure gradient (explicit) 
    116                               CALL iom_put( "vtrd_spgexp", pvtrd ) 
    117       CASE( jpdyn_spgflt );   CALL iom_put( "utrd_spgflt", putrd ) ! surface pressure gradient (filtered) 
    118                               CALL iom_put( "vtrd_spgflt", pvtrd ) 
    119115      CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity 
    120116                              CALL iom_put( "vtrd_pvo", pvtrd ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

    r5901 r6079  
    120120         CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg", zke )    ! hydrostatic pressure gradient 
    121121         CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg", zke )    ! surface pressure gradient 
    122          CASE( jpdyn_spgexp );   CALL iom_put( "ketrd_spgexp", zke ) ! surface pressure gradient (explicit) 
    123          CASE( jpdyn_spgflt );   CALL iom_put( "ketrd_spgflt", zke ) ! surface pressure gradient (filter) 
    124122         CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo", zke )    ! planetary vorticity 
    125123         CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo", zke )    ! relative  vorticity     (or metric term) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r5901 r6079  
    2121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s] 
    2222   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s] 
    23    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ua_sv,  va_sv          !: Saved trends (time spliting) [m/s2] 
    2423   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s] 
    2524   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1] 
     
    3837   ! 
    3938   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient 
     39 
     40   !! Specific to split explicit free surface (allocated in dynspg_ts module): 
     41   ! 
     42   !! Time filtered arrays at baroclinic time step: 
     43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv    ! Advection vel. at "now" barocl. step 
     44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b , vb2_b     ! Half step fluxes (ln_bt_fw=T) 
     45#if defined key_agrif 
     46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_i_b,  vb2_i_b ! Half step time integrated fluxes  
     47#endif 
     48 
     49   !! Arrays at barotropic time step:                   ! bef before !   before   !    now     !   after    ! 
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ubb_e    ,    ub_e    ,    un_e    ,    ua_e  
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vbb_e    ,    vb_e    ,    vn_e    ,    va_e  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshbb_e  ,    sshb_e  ,    sshn_e  ,    ssha_e  
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hu_e 
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hv_e                                                 
     55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hur_e 
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::                              hvr_e  
    4057 
    4158   !! interpolated gradient (only used in zps case) 
     
    85102      ! 
    86103      ALLOCATE( ub   (jpi,jpj,jpk)      , un   (jpi,jpj,jpk)      , ua(jpi,jpj,jpk)       ,     & 
    87          &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &       
    88          &      ua_sv(jpi,jpj,jpk)      , va_sv(jpi,jpj,jpk)      ,                             &       
     104         &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &           
    89105         &      wn   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             & 
    90106         &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     & 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5901 r6079  
    2828   !!            3.7  !  2014-10  (G. Madec)  LDF simplication  
    2929   !!             -   !  2014-12  (G. Madec) remove KPP scheme 
     30   !!             -   !  2015-11  (J. Chanut) free surface simplification 
    3031   !!---------------------------------------------------------------------- 
    3132 
     
    179180 
    180181      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    181       !  Ocean dynamics : hdiv, ssh, e3, wn 
    182       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     182      !  Ocean dynamics : hdiv, ssh, e3, u, v, w 
     183      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     184 
    183185                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor) 
    184186      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    185187                         CALL wzv           ( kstp )  ! now cross-level velocity  
    186188 
    187       IF( lk_dynspg_ts ) THEN  
    188           ! In case the time splitting case, update almost all momentum trends here: 
    189           ! Note that the computation of vertical velocity above, hence "after" sea level 
    190           ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    191189!!gm : why also here ???? 
    192             IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations 
     190      IF(ln_sto_eos  )  CALL sto_pts( tsn )                             ! Random T/S fluctuations 
    193191!!gm 
    194                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    195                              
     192                         CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
     193 
     194!!jc: fs simplification 
     195!!jc: lines below are useless if lk_vvl=T. Keep them here (which maintains a bug if lk_vvl=F and ln_zps=T, cf ticket #1636)  
     196!!                                         but ensures reproductible results 
     197!!                                         with previous versions using split-explicit free surface           
    196198            IF( ln_zps .AND. .NOT. ln_isfcav)   &                           ! Partial steps: bottom before horizontal gradient 
    197199               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &  ! of t, s, rd at the last ocean level 
     
    201203               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    202204               &                                               grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) 
    203  
    204                                   ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
    205                                   va(:,:,:) = 0._wp 
    206           IF(  lk_asminc .AND. ln_asmiau .AND. & 
    207              & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment 
    208           IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends 
    209                                   CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
    210                                   CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
    211                                   CALL dyn_ldf      ( kstp )   ! lateral mixing 
    212 #if defined key_agrif 
    213           IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge 
    214 #endif 
    215                                   CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    216                                   CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    217  
    218                                   ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated at this stage) 
    219                                   va_sv(:,:,:) = va(:,:,:) 
    220  
    221                                   CALL div_hor( kstp )         ! Horizontal divergence  (2nd call in time-split case) 
    222           IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    223                                   CALL wzv           ( kstp )  ! now cross-level velocity  
    224       ENDIF 
    225  
    226       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    227       ! diagnostics and outputs             (ua, va, tsa used as workspace) 
    228       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    229       IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats 
    230       IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth) 
    231       IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics 
    232       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports 
    233       IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag 
    234       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    235                          CALL dia_wri( kstp )         ! ocean model: outputs 
    236       ! 
    237       IF( ln_crs     )   CALL crs_fld( kstp )         ! ocean model: online field coarsening & output 
     205!!jc: fs simplification 
     206                             
     207                         ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
     208                         va(:,:,:) = 0._wp 
     209 
     210      IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
     211                         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment 
     212      IF( lk_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends 
     213#if defined key_agrif 
     214      IF(.NOT. Agrif_Root())  &  
     215               &         CALL Agrif_Sponge_dyn        ! momentum sponge 
     216#endif 
     217                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
     218                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
     219                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
     220                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
     221                         CALL dyn_spg       ( kstp )  ! surface pressure gradient 
     222 
     223                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well 
     224      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
     225                            CALL div_hor    ( kstp )              ! Horizontal divergence  (2nd call in time-split case) 
     226         IF( lk_vvl )       CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
     227                            CALL wzv        ( kstp )              ! now cross-level velocity  
     228      ENDIF 
     229 
     230                         CALL dyn_bfr       ( kstp )  ! bottom friction 
     231                         CALL dyn_zdf       ( kstp )  ! vertical diffusion 
     232 
     233      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     234      ! diagnostics and outputs              
     235      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     236      IF( lk_floats  )   CALL flo_stp       ( kstp )  ! drifting Floats 
     237      IF( lk_diahth  )   CALL dia_hth       ( kstp )  ! Thermocline depth (20 degres isotherm depth) 
     238      IF(.NOT.ln_cpl )   CALL dia_fwb       ( kstp )  ! Fresh water budget diagnostics 
     239      IF( lk_diadct  )   CALL dia_dct       ( kstp )  ! Transports 
     240      IF( lk_diaar5  )   CALL dia_ar5       ( kstp )  ! ar5 diag 
     241      IF( lk_diaharm )   CALL dia_harm      ( kstp )  ! Tidal harmonic analysis 
     242                         CALL dia_wri       ( kstp )  ! ocean model: outputs 
     243      ! 
     244      IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output 
    238245 
    239246#if defined key_top 
     
    241248      ! Passive Tracer Model 
    242249      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    243                          CALL trc_stp( kstp )         ! time-stepping 
    244 #endif 
    245  
    246       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    247       ! Active tracers                              (ua, va used as workspace) 
    248       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    249                              tsa(:,:,:,:) = 0._wp           ! set tracer trends to zero 
     250                         CALL trc_stp       ( kstp )  ! time-stepping 
     251#endif 
     252 
     253      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     254      ! Active tracers                               
     255      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     256                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
    250257 
    251258      IF(  lk_asminc .AND. ln_asmiau .AND. & 
    252          & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment 
    253                              CALL tra_sbc    ( kstp )       ! surface boundary condition 
    254       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
    255       IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux 
    256       IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
    257       IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends 
    258       IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
    259                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection 
    260                              CALL tra_ldf    ( kstp )       ! lateral mixing 
     259         & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
     260                         CALL tra_sbc       ( kstp )  ! surface boundary condition 
     261      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr 
     262      IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux 
     263      IF( lk_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme 
     264      IF( ln_tradmp  )   CALL tra_dmp       ( kstp )  ! internal damping trends 
     265      IF( lk_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends 
     266#if defined key_agrif 
     267      IF(.NOT. Agrif_Root())  &  
     268               &         CALL Agrif_Sponge_tra        ! tracers sponge 
     269#endif 
     270                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection 
     271                         CALL tra_ldf       ( kstp )  ! lateral mixing 
    261272 
    262273!!gm : why CALL to dia_ptr has been moved here??? (use trends info?) 
    263       IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics 
     274      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics 
    264275!!gm 
    265  
    266 #if defined key_agrif 
    267       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge 
    268 #endif 
    269                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields 
    270  
    271       IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos) 
    272          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    273                              CALL tra_nxt( kstp )                ! tracer fields at next time step 
    274 !!gm : why again a call to sto_pts ??? 
    275             IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations 
    276 !!gm 
    277                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    278             IF( ln_zps .AND. .NOT. ln_isfcav)                                & 
    279                &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    280                &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    281             IF( ln_zps .AND.       ln_isfcav)                                & 
    282                &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv, gtui, gtvi,  &    ! Partial steps for top/bottom cells 
    283                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    284                &                                                grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) 
    285       ELSE                                                  ! centered hpg  (eos then time stepping) 
    286          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    287 !!gm : why again a call to sto_pts ??? 
    288             IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations 
    289 !!gm 
    290                              CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
    291          IF( ln_zps .AND. .NOT. ln_isfcav)                                   & 
    292                &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: bottom before horizontal gradient 
    293                &                                           rhd, gru , grv    )    ! of t, s, rd at the last ocean level 
    294          IF( ln_zps .AND.       ln_isfcav)                                   &  
    295                &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,   &    ! Partial steps for top/bottom cells 
    296                &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    297                &                                    grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    298          ENDIF 
    299          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    300                              CALL tra_nxt( kstp )                ! tracer fields at next time step 
    301       ENDIF 
    302  
    303       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    304       ! Dynamics                                    (tsa used as workspace) 
    305       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    306       IF( lk_dynspg_ts   )  THEN 
    307                                                              ! revert to previously computed momentum tendencies 
    308                                                              ! (not using ua, va as temporary arrays during tracers' update could avoid that) 
    309                                ua(:,:,:) = ua_sv(:,:,:) 
    310                                va(:,:,:) = va_sv(:,:,:) 
    311  
    312                                CALL dyn_bfr( kstp )         ! bottom friction 
    313                                CALL dyn_zdf( kstp )         ! vertical diffusion 
    314       ELSE 
    315                                ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
    316                                va(:,:,:) = 0._wp 
    317  
    318         IF(  lk_asminc .AND. ln_asmiau .AND. & 
    319            & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
    320         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
    321         IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
    322                                CALL dyn_adv( kstp )         ! advection (vector or flux form) 
    323                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis 
    324                                CALL dyn_ldf( kstp )         ! lateral mixing 
    325 #if defined key_agrif 
    326         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge 
    327 #endif 
    328                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    329                                CALL dyn_bfr( kstp )         ! bottom friction 
    330                                CALL dyn_zdf( kstp )         ! vertical diffusion 
    331                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    332       ENDIF 
    333                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
    334  
    335                                CALL ssh_swp( kstp )         ! swap of sea surface height 
    336       IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
     276                         CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields 
     277      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection 
     278 
     279      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     280      ! Set boundary conditions and Swap 
     281      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     282!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap  
     283!!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.  
     284!!    If so:  
     285!!    (i) no need to call agrif update at initialization time 
     286!!    (ii) no need to update "before" fields  
     287!! 
     288!!    Apart from creating new tra_swp/dyn_swp routines, this however:  
     289!!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between  
     290!!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",  
     291!!    e.g. a shift of the feedback interface inside child domain.  
     292!!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 
     293!!    place. 
     294!!  
     295!!jc2: dynnxt must be the latest call. fse3t_b are indeed updated in that routine 
     296                         CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap 
     297                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap 
     298                         CALL ssh_swp       ( kstp )  ! swap of sea surface height 
     299      IF( lk_vvl     )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    337300      ! 
    338301 
    339302!!gm : This does not only concern the dynamics ==>>> add a new title 
    340303!!gm2: why ouput restart before AGRIF update? 
    341       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
     304!! 
     305!!jc: That would be better, but see comment above 
     306!! 
     307      IF( lrst_oce   )   CALL rst_write     ( kstp )  ! write output ocean restart file 
    342308 
    343309#if defined key_agrif 
     
    345311      ! AGRIF 
    346312      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    347                                CALL Agrif_Integrate_ChildGrids( stp )   
    348  
    349       IF ( Agrif_NbStepint().EQ.0 ) THEN 
    350                                CALL Agrif_Update_Tra()      ! Update active tracers 
    351                                CALL Agrif_Update_Dyn()      ! Update momentum 
    352       ENDIF 
    353 #endif 
    354       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics 
    355       IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     313                         CALL Agrif_Integrate_ChildGrids( stp )   
     314 
     315      IF ( Agrif_NbStepint().EQ.0 ) THEN              ! AGRIF Update  
     316!!jc in fact update i useless at last time step, but do it for global diagnostics 
     317                         CALL Agrif_Update_Tra()      ! Update active tracers 
     318                         CALL Agrif_Update_Dyn()      ! Update momentum 
     319      ENDIF 
     320#endif 
     321      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics 
     322      IF( lk_diaobs  )   CALL dia_obs       ( kstp )  ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    356323 
    357324      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    358325      ! Control 
    359326      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    360                                CALL stp_ctl( kstp, indic ) 
    361       IF( indic < 0        )   THEN 
    362                                CALL ctl_stop( 'step: indic < 0' ) 
    363                                CALL dia_wri_state( 'output.abort', kstp ) 
    364       ENDIF 
    365       IF( kstp == nit000   )   THEN 
     327                         CALL stp_ctl       ( kstp, indic ) 
     328      IF( indic < 0  )   THEN 
     329                         CALL ctl_stop( 'step: indic < 0' ) 
     330                         CALL dia_wri_state( 'output.abort', kstp ) 
     331      ENDIF 
     332      IF( kstp == nit000 )   THEN 
    366333                 CALL iom_close( numror )     ! close input  ocean restart file 
    367334         IF(lwm) CALL FLUSH    ( numond )     ! flush output namelist oce 
     
    374341      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    375342!!gm why lk_oasis and not lk_cpl ???? 
    376       IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges 
     343      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges 
    377344      ! 
    378345#if defined key_iomput 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r5901 r6079  
    4040   USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine) 
    4141   USE dynzdf           ! vertical diffusion               (dyn_zdf routine) 
    42    USE dynspg_oce       ! surface pressure gradient        (dyn_spg routine) 
    4342   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    4443 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/stpctl.F90

    r3294 r6079  
    1616   USE oce             ! ocean dynamics and tracers variables 
    1717   USE dom_oce         ! ocean space and time domain variables  
    18    USE sol_oce         ! ocean space and time domain variables  
    1918   USE in_out_manager  ! I/O manager 
    2019   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2120   USE lib_mpp         ! distributed memory computing 
    22    USE dynspg_oce      ! pressure gradient schemes  
    2321   USE c1d             ! 1D vertical configuration 
    2422 
     
    4341      !! ** Method  : - Save the time step in numstp 
    4442      !!              - Print it each 50 time steps 
    45       !!              - Print solver statistics in numsol  
    46       !!              - Stop the run IF problem for the solver ( indec < 0 ) 
     43      !!              - Stop the run IF problem ( indic < 0 ) 
    4744      !! 
    4845      !! ** Actions :   'time.step' file containing the last ocean time-step 
     
    5047      !!---------------------------------------------------------------------- 
    5148      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
    52       INTEGER, INTENT( inout ) ::   kindic  ! indicator of solver convergence 
     49      INTEGER, INTENT( inout ) ::   kindic  ! error indicator 
    5350      !! 
    5451      INTEGER  ::   ji, jj, jk              ! dummy loop indices 
     
    143140      IF( lk_c1d )  RETURN          ! No log file in case of 1D vertical configuration 
    144141 
    145       ! log file (solver or ssh statistics) 
    146       ! -------- 
    147       IF( lk_dynspg_flt ) THEN      ! elliptic solver statistics (if required) 
    148          ! 
    149          IF(lwp) WRITE(numsol,9200) kt, niter, res, SQRT(epsr)/eps       ! Solver 
    150          ! 
    151          IF( kindic < 0 .AND. zsmin > 0.e0 .AND. zumax <= 20.e0 ) THEN   ! create a abort file if problem found  
    152             IF(lwp) THEN 
    153                WRITE(numout,*) ' stpctl: the elliptic solver DO not converge or explode' 
    154                WRITE(numout,*) ' ====== ' 
    155                WRITE(numout,9200) kt, niter, res, sqrt(epsr)/eps 
    156                WRITE(numout,*) 
    157                WRITE(numout,*) ' stpctl: output of last fields' 
    158                WRITE(numout,*) ' ======  ' 
    159             ENDIF 
    160          ENDIF 
    161          ! 
    162       ELSE                                   !* ssh statistics (and others...) 
    163          IF( kt == nit000 .AND. lwp ) THEN   ! open ssh statistics file (put in solver.stat file) 
    164             CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    165          ENDIF 
    166          ! 
    167          zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 
    168          IF( lk_mpp )   CALL mpp_sum( zssh2 )      ! sum over the global domain 
    169          ! 
    170          IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin      ! ssh statistics 
    171          ! 
     142      ! log file (ssh statistics) 
     143      ! --------                                   !* ssh statistics (and others...) 
     144      IF( kt == nit000 .AND. lwp ) THEN   ! open ssh statistics file (put in solver.stat file) 
     145         CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    172146      ENDIF 
     147      ! 
     148      zssh2 = SUM( sshn(:,:) * sshn(:,:) * tmask_i(:,:) ) 
     149      IF( lk_mpp )   CALL mpp_sum( zssh2 )      ! sum over the global domain 
     150      ! 
     151      IF(lwp) WRITE(numsol,9300) kt, zssh2, zumax, zsmin      ! ssh statistics 
     152      ! 
    173153 
    1741549200  FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC/diawri.F90

    r5901 r6079  
    2626   USE dom_oce         ! ocean space and time domain 
    2727   USE zdf_oce         ! ocean vertical physics 
    28    USE sol_oce         ! solver variables 
    2928   USE sbc_oce         ! Surface boundary condition: ocean fields 
    3029   USE sbc_ice         ! Surface boundary condition: ice fields 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/TOP_SRC

    • Property svn:mergeinfo deleted
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/TOP_SRC/trcsub.F90

    r5901 r6079  
    502502      z1_rau0 = 0.5 / rau0 
    503503      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    504 #if ! defined key_dynspg_ts 
     504 
    505505      ! These lines are not necessary with time splitting since 
    506506      ! boundary condition on sea level is set during ts loop 
     
    512512      CALL lbc_lnk( ssha, 'T', 1. )  
    513513#endif 
    514 #endif 
    515  
    516514 
    517515      !                                           !------------------------------! 
Note: See TracChangeset for help on using the changeset viewer.