Changeset 12614


Ignore:
Timestamp:
2020-03-26T15:59:52+01:00 (7 months ago)
Author:
gm
Message:

first Shallow Water Eq. update

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater
Files:
22 added
13 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90

    r12529 r12614  
    896896         IF ( kt == nitdin_r ) THEN 
    897897            ! 
    898             l_1st_euler = 0              ! Force Euler forward step 
     898            l_1st_euler = .TRUE.              ! Force Euler forward step 
    899899            ! 
    900900            ! Sea-ice : SI3 case 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90

    r12529 r12614  
    119119      REAL(wp)::   zztmp , zztmpx   ! local scalar 
    120120      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     121      REAL(wp)::   ze3              !   -      - 
    121122      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
    122123      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
     
    146147         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
    147148      ELSE 
    148          CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    149       ENDIF 
     149         CALL iom_put( "ssh" , ssh(:,:,Kmm) )                          ! sea surface height 
     150      ENDIF 
     151 
     152!!an  
     153        IF( iom_use("ht") )   &                  ! water column at t-point 
     154         CALL iom_put( "ht" , ht_0(:,:) + ssh(:,:,Kmm) ) 
     155      ! 
     156      IF( iom_use("hu") )   THEN                   ! water column at u-point 
     157         z2d(:,:) = 0._wp 
     158         DO_2D_10_10   
     159            z2d(ji,jj) = 0.5_wp * ( e3t(ji  ,jj,1,Kmm) * e1e2t(ji  ,jj)  & 
     160               &                  + e3t(ji+1,jj,1,Kmm) * e1e2t(ji+1,jj)  ) * r1_e1e2u(ji,jj) 
     161         END_2D 
     162         CALL lbc_lnk( 'diawri', z2d, 'U', 1._wp ) 
     163         CALL iom_put( "hu", z2d )  
     164      ENDIF 
     165     !                     
     166      IF( iom_use("hv") )   THEN                  ! water column at v-point 
     167         z2d(:,:) = 0._wp 
     168         DO_2D_10_10   
     169            z2d(ji,jj) = 0.5_wp * (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1)   & 
     170              &                    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  )   ) * r1_e1e2v(ji,jj) 
     171         END_2D 
     172         CALL lbc_lnk( 'diawri', z2d, 'V', 1._wp ) 
     173         CALL iom_put( "hv", z2d )     
     174      ENDIF              
     175      ! 
     176      IF( iom_use("hf") )  THEN                    ! water column at f-point 
     177         z2d(:,:) = 0._wp 
     178         DO_2D_10_10   
     179            z2d(ji,jj) = 0.25_wp * (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     180              &                     + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     181         END_2D 
     182         CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp ) 
     183         CALL iom_put( "hf", z2d )    
     184      ENDIF               
     185!!an 
    150186 
    151187      IF( iom_use("wetdep") )   &                  ! wet depth 
    152188         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    153        
    154       CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
    155       CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    156       IF ( iom_use("sbt") ) THEN 
    157          DO_2D_11_11 
    158             ikbot = mbkt(ji,jj) 
    159             z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
    160          END_2D 
    161          CALL iom_put( "sbt", z2d )                ! bottom temperature 
    162       ENDIF 
    163        
    164       CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity 
    165       CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    166       IF ( iom_use("sbs") ) THEN 
    167          DO_2D_11_11 
    168             ikbot = mbkt(ji,jj) 
    169             z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
    170          END_2D 
    171          CALL iom_put( "sbs", z2d )                ! bottom salinity 
    172       ENDIF 
    173189 
    174190      IF ( iom_use("taubot") ) THEN                ! bottom stress 
     
    186202         CALL iom_put( "taubot", z2d )            
    187203      ENDIF 
    188           
    189       CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current 
    190       CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    191       IF ( iom_use("sbu") ) THEN 
    192          DO_2D_11_11 
    193             ikbot = mbku(ji,jj) 
    194             z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     204      !   
     205      CALL iom_put(  "ssu" , uu(:,:,1,Kmm) )            ! surface i-current 
     206      CALL iom_put(  "ssv" , vv(:,:,1,Kmm) )            ! surface j-current 
     207      CALL iom_put(  "woce", ww )                       ! vertical velocity 
     208      ! 
     209       
     210      IF ( iom_use("sKE") ) THEN                        ! surface kinetic energy at T point 
     211         z2d(:,:) = 0._wp 
     212         DO_2D_00_00 
     213            z2d(ji,jj) =    0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  & 
     214               &                      + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm)  & 
     215               &                      + vv(ji,jj  ,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,1,Kmm)  &  
     216               &                      + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm)  )  & 
     217               &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * tmask(ji,jj,1) 
    195218         END_2D 
    196          CALL iom_put( "sbu", z2d )                ! bottom i-current 
     219         ! 
     220         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
     221         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )    
     222                            
     223      ENDIF 
     224            
     225!!an sKEf 
     226      IF ( iom_use("sKEf") ) THEN                        ! surface kinetic energy at F point 
     227         z2d(:,:) = 0._wp                                ! CAUTION : only valid in SWE, not with bathymetry 
     228         DO_2D_00_00 
     229            z2d(ji,jj) =    0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  & 
     230               &                      + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm)  & 
     231               &                      + vv(ji  ,jj,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji  ,jj) * e3v(ji  ,jj,1,Kmm)  &  
     232               &                      + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm)  )  & 
     233               &                    * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj) 
     234         END_2D 
     235         ! 
     236         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     237         CALL iom_put( "sKEf", z2d )                      
     238      ENDIF 
     239!!an 
     240      ! 
     241      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
     242      ! 
     243      ! Output of vorticity terms 
     244      IF ( iom_use("relvor")    .OR. iom_use("plavor")    .OR.   & 
     245         & iom_use("relpotvor") .OR. iom_use("abspotvor") .OR.   & 
     246         & iom_use("Ens")                                        ) THEN 
     247         ! 
     248         z2d(:,:) = 0._wp  
     249         ze3 = 0._wp  
     250         DO_2D_10_10 
     251            z2d(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    & 
     252            &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm)  ) * r1_e1e2f(ji,jj) 
     253         END_2D 
     254         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     255         CALL iom_put( "relvor", z2d )                  ! relative vorticity ( zeta )  
     256         ! 
     257         CALL iom_put( "plavor", ff_f )                 ! planetary vorticity ( f ) 
     258         ! 
     259         DO_2D_10_10   
     260            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     261              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     262            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     263            ELSE                      ;   ze3 = 0._wp 
     264            ENDIF 
     265            z2d(ji,jj) = ze3 * z2d(ji,jj)  
     266         END_2D 
     267         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     268         CALL iom_put( "relpotvor", z2d )                  ! relative potential vorticity (zeta/h) 
     269         ! 
     270         DO_2D_10_10 
     271            ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    & 
     272              &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj) 
     273            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     274            ELSE                      ;   ze3 = 0._wp 
     275            ENDIF 
     276            z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)  
     277         END_2D 
     278         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     279         CALL iom_put( "abspotvor", z2d )                  ! absolute potential vorticity ( q ) 
     280         ! 
     281         DO_2D_10_10   
     282            z2d(ji,jj) = 0.5_wp * z2d(ji,jj)  * z2d(ji,jj)  
     283         END_2D 
     284         CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 
     285         CALL iom_put( "Ens", z2d )                        ! potential enstrophy ( 1/2*q2 ) 
     286         ! 
    197287      ENDIF 
    198288       
    199       CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current 
    200       CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    201       IF ( iom_use("sbv") ) THEN 
    202          DO_2D_11_11 
    203             ikbot = mbkv(ji,jj) 
    204             z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
    205          END_2D 
    206          CALL iom_put( "sbv", z2d )                ! bottom j-current 
    207       ENDIF 
    208  
    209       IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
    210       ! 
    211       CALL iom_put( "woce", ww )                   ! vertical velocity 
    212       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    213          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    214          z2d(:,:) = rho0 * e1e2t(:,:) 
    215          DO jk = 1, jpk 
    216             z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:) 
    217          END DO 
    218          CALL iom_put( "w_masstr" , z3d )   
    219          IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    220       ENDIF 
    221       ! 
    222       IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    223  
    224       CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
    225       CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef. 
    226       CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef. 
    227  
    228       IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) ) 
    229       IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    230  
    231       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    232          DO_2D_00_00 
    233             zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
    234             zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 
    235             zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 
    236             z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    237                &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    238          END_2D 
    239          CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    240          CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
    241          z2d(:,:) = SQRT( z2d(:,:) ) 
    242          CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
    243       ENDIF 
    244           
    245       ! heat and salt contents 
    246       IF( iom_use("heatc") ) THEN 
    247          z2d(:,:)  = 0._wp  
    248          DO_3D_11_11( 1, jpkm1 ) 
    249             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    250          END_3D 
    251          CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    252       ENDIF 
    253  
    254       IF( iom_use("saltc") ) THEN 
    255          z2d(:,:)  = 0._wp  
    256          DO_3D_11_11( 1, jpkm1 ) 
    257             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    258          END_3D 
    259          CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    260       ENDIF 
    261       ! 
    262       IF ( iom_use("eken") ) THEN 
    263          z3d(:,:,jpk) = 0._wp  
    264          DO_3D_00_00( 1, jpkm1 ) 
    265             zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    266             z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
    267                &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
    268                &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
    269                &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    270          END_3D 
    271          CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    272          CALL iom_put( "eken", z3d )                 ! kinetic energy 
    273       ENDIF 
    274       ! 
    275       CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
    276       ! 
    277       IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
    278          z3d(:,:,jpk) = 0.e0 
    279          z2d(:,:) = 0.e0 
    280          DO jk = 1, jpkm1 
    281             z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk) 
    282             z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    283          END DO 
    284          CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction 
    285          CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum 
    286       ENDIF 
    287        
    288       IF( iom_use("u_heattr") ) THEN 
    289          z2d(:,:) = 0._wp  
    290          DO_3D_00_00( 1, jpkm1 ) 
    291             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
    292          END_3D 
    293          CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    294          CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
    295       ENDIF 
    296  
    297       IF( iom_use("u_salttr") ) THEN 
    298          z2d(:,:) = 0.e0  
    299          DO_3D_00_00( 1, jpkm1 ) 
    300             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
    301          END_3D 
    302          CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    303          CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
    304       ENDIF 
    305  
    306        
    307       IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN 
    308          z3d(:,:,jpk) = 0.e0 
    309          DO jk = 1, jpkm1 
    310             z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    311          END DO 
    312          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
    313       ENDIF 
    314        
    315       IF( iom_use("v_heattr") ) THEN 
    316          z2d(:,:) = 0.e0  
    317          DO_3D_00_00( 1, jpkm1 ) 
    318             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
    319          END_3D 
    320          CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    321          CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
    322       ENDIF 
    323  
    324       IF( iom_use("v_salttr") ) THEN 
    325          z2d(:,:) = 0._wp  
    326          DO_3D_00_00( 1, jpkm1 ) 
    327             z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
    328          END_3D 
    329          CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    330          CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
    331       ENDIF 
    332  
    333       IF( iom_use("tosmint") ) THEN 
    334          z2d(:,:) = 0._wp 
    335          DO_3D_00_00( 1, jpkm1 ) 
    336             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    337          END_3D 
    338          CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    339          CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
    340       ENDIF 
    341       IF( iom_use("somint") ) THEN 
    342          z2d(:,:)=0._wp 
    343          DO_3D_00_00( 1, jpkm1 ) 
    344             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    345          END_3D 
    346          CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    347          CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
    348       ENDIF 
    349  
    350       CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
    351       ! 
    352        
    353       IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
    354  
     289      ! 
    355290      IF( ln_timing )   CALL timing_stop('dia_wri') 
    356291      ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90

    r12529 r12614  
    151151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  !: u-depth              [m] 
    152152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  !: v-depth              [m] 
     153   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0  !: f-depth              [m] 
     154   !                                                      !  inverse of reference heights of water column 
     155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_ht_0  !: t-depth              [1/m] 
     156   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hu_0  !: u-depth              [1/m] 
     157   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hv_0  !: v-depth              [1/m] 
     158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   r1_hf_0  !: f-depth              [1/m] 
     159    
    153160                                                          ! time-dependent heights of water column 
    154161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht                     !: height of water column at T-points [m] 
     
    178185   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF) 
    179186 
    180    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
     187   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts 
    181188   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    182189   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
     
    268275         &      e3uw_0(jpi,jpj,jpk)     , e3vw_0(jpi,jpj,jpk)     ,         & 
    269276         &      e3uw  (jpi,jpj,jpk,jpt) , e3vw  (jpi,jpj,jpk,jpt) ,    STAT=ierr(5) )                        
    270          ! 
    271       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj)    , hv_0(jpi,jpj)     ,                                             & 
    272          &      ht  (jpi,jpj) , hu(  jpi,jpj,jpt), hv(  jpi,jpj,jpt) , r1_hu(jpi,jpj,jpt) , r1_hv(jpi,jpj,jpt) ,   & 
    273          &                      STAT=ierr(6)  ) 
     277         !     
     278      ALLOCATE(    ht_0(jpi,jpj) ,    hu_0(jpi,jpj)     ,    hv_0(jpi,jpj)     ,    hf_0(jpi,jpj) ,     & 
     279         &      r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj)     , r1_hv_0(jpi,jpj)     , r1_hf_0(jpi,jpj) ,     & 
     280         &         ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt) ,    hv  (jpi,jpj,jpt)                    ,     & 
     281         &                         r1_hu  (jpi,jpj,jpt) , r1_hv  (jpi,jpj,jpt)                    , STAT=ierr(6)  ) 
    274282         ! 
    275283      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  )  
     
    277285      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) ) 
    278286         ! 
    279       ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
    280          &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     & 
    281          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
     287      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                           &  
     288         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) ,     & 
     289         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj)                    , STAT=ierr(9) ) 
    282290         ! 
    283291      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90

    r12529 r12614  
    143143 
    144144      CALL dom_msk( ik_top, ik_bot )    ! Masks 
    145       ! 
    146       ht_0(:,:) = 0._wp  ! Reference ocean thickness 
     145 
     146      ht_0(:,:) = 0._wp                 ! Reference ocean thickness 
    147147      hu_0(:,:) = 0._wp 
    148148      hv_0(:,:) = 0._wp 
    149       DO jk = 1, jpk 
     149      hf_0(:,:) = 0._wp 
     150      DO jk = 1, jpkm1 
    150151         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    151152         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 
    152153         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 
     154         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * ssfmask(:,:)     ! CAUTION : only valid in SWE, not with bathymetry 
    153155      END DO 
    154       ! 
     156      !                                 ! Inverse of reference ocean thickness 
     157      r1_ht_0(:,:) =  ssmask(:,:) / ( ht_0(:,:) + 1._wp -  ssmask(:,:) ) 
     158      r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) 
     159      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
     160      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) ) 
     161       
    155162      !           !==  time varying part of coordinate system  ==! 
    156163      ! 
     
    190197      ! 
    191198      IF( ln_meshmask    )   CALL dom_wri       ! Create a domain file 
     199 
    192200      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control 
    193201      ! 
     
    200208         WRITE(numout,*)  
    201209      ENDIF 
     210       
    202211      ! 
    203212   END SUBROUTINE dom_init 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dommsk.F90

    r12529 r12614  
    192192      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 
    193193      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 
    194  
     194!!an  
     195      ! ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 
     196      DO_2D_10_10 
     197         ssfmask(ji,jj) = MAX(  tmask(ji,jj+1,1), tmask(ji+1,jj+1,1),  &  
     198            &                   tmask(ji,jj  ,1), tmask(ji+1,jj  ,1)   ) 
     199      END_2D 
     200      CALL lbc_lnk( 'dommsk', ssfmask, 'F', 1._wp )    
     201 
     202!!an 
    195203 
    196204      ! Interior domain mask  (used for global sum) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90

    r12529 r12614  
    166166      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    167167      !                                ! Horizontal interpolation of e3t 
    168       CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
    169       CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    170       CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
    171       CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    172       CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
     168      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
     169      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     170      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     171      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     172      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    173173      !                                ! Vertical interpolation of e3t,u,v  
    174       CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
    175       CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
    176       CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW 
    177       CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
    178       CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW 
    179       CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
     174      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
     175      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
     176      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW 
     177      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     178      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW 
     179      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    180180 
    181181      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 
     
    549549      ! *********************************** ! 
    550550 
    551       CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 
    552       CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 
     551      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 
     552      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 
    553553 
    554554      ! *********************************** ! 
     
    633633      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    634634       
    635       CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
     635      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F'  ) 
    636636       
    637637      ! Vertical scale factor interpolations 
    638       CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
    639       CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
    640       CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
    641       CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  ) 
    642       CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
    643       CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
     638      CALL dom_vvl_interpol( ssh(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
     639      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     640      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     641      CALL dom_vvl_interpol( ssh(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  ) 
     642      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     643      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    644644 
    645645      ! t- and w- points depth (set the isf depth as it is in the initial step) 
     
    674674 
    675675 
    676    SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
     676   SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp ) 
    677677      !!--------------------------------------------------------------------- 
    678678      !!                  ***  ROUTINE dom_vvl__interpol  *** 
     
    684684      !!                - vertical interpolation: simple averaging 
    685685      !!---------------------------------------------------------------------- 
    686       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated 
    687       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3 
    688       CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors 
    689       !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    690       ! 
    691       INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
    692       REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F 
    693       !!---------------------------------------------------------------------- 
    694       ! 
    695       IF(ln_wd_il) THEN 
    696         zlnwd = 1.0_wp 
    697       ELSE 
    698         zlnwd = 0.0_wp 
    699       END IF 
    700       ! 
    701       SELECT CASE ( pout )    !==  type of interpolation  ==! 
     686      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pssh    ! input e3   NOT used here (ssh is used instead) 
     687      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m] 
     688      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' ) 
     689      ! 
     690      INTEGER ::   ji, jj, jk                 ! dummy loop indices 
     691      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace 
     692      !!---------------------------------------------------------------------- 
     693      ! 
     694      SELECT CASE ( cdp )     !==  type of interpolation  ==! 
    702695         ! 
    703696      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    704          DO_3D_10_10( 1, jpk ) 
    705             pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    706                &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    707                &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    708          END_3D 
    709          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    710          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
     697         DO_2D_00_00 
     698            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  & 
     699               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj) 
     700         END_2D 
     701         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp ) 
     702         ! 
     703         DO jk = 1, jpkm1 
     704            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 
     705         END DO 
    711706         ! 
    712707      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    713          DO_3D_10_10( 1, jpk ) 
    714             pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    715                &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    716                &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    717          END_3D 
    718          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    719          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
     708         DO_2D_00_00 
     709            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  & 
     710               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 
     711         END_2D 
     712         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp ) 
     713         ! 
     714         DO jk = 1, jpkm1 
     715            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 
     716         END DO 
    720717         ! 
    721718      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    722          DO_3D_10_10( 1, jpk ) 
    723             pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    724                &                       *    r1_e1e2f(ji,jj)                                                  & 
    725                &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    726                &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    727          END_3D 
    728          CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    729          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     719         DO_2D_10_10 
     720            zc3(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  & 
     721               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  & 
     722               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  & 
     723               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 
     724         END_2D 
     725         CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 
     726         ! 
     727         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh 
     728            e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + zc3(:,:) ) 
     729         END DO 
    730730         ! 
    731731      CASE( 'W' )                   !* from T- to W-point : vertical simple mean 
    732          ! 
    733          pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 
    734          ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 
    735 !!gm BUG? use here wmask in case of ISF ?  to be checked 
    736          DO jk = 2, jpk 
    737             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   & 
    738                &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               & 
    739                &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     & 
    740                &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    741          END DO 
    742          ! 
    743       CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean 
    744          ! 
    745          pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 
    746          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    747 !!gm BUG? use here wumask in case of ISF ?  to be checked 
    748          DO jk = 2, jpk 
    749             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
    750                &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              & 
    751                &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    752                &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    753          END DO 
    754          ! 
    755       CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean 
    756          ! 
    757          pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 
    758          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    759 !!gm BUG? use here wvmask in case of ISF ?  to be checked 
    760          DO jk = 2, jpk 
    761             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
    762                &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              & 
    763                &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    764                &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    765          END DO 
     732         zc3(:,:) = pssh(:,:) * r1_ht_0(:,:) 
     733         ! 
     734         DO jk = 1, jpk 
     735            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 
     736         END DO 
     737         ! 
     738      CASE( 'UW' )                  !* from U- to UW-point 
     739         ! 
     740         DO_2D_00_00 
     741            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  & 
     742               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj) 
     743         END_2D 
     744         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp ) 
     745         ! 
     746         DO jk = 1, jpk 
     747            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 
     748         END DO 
     749      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean 
     750         ! 
     751         DO_2D_00_00 
     752            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  & 
     753               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 
     754         END_2D 
     755         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp ) 
     756          ! 
     757         DO jk = 1, jpk 
     758            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 
     759         END DO 
     760         ! 
    766761      END SELECT 
    767762      ! 
     
    878873                  ! Wetting and drying test case 
    879874                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    880                   ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     875!!an                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    881876                  ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
    882877                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     
    923918!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    924919                ssh(:,:,Kmm)=0._wp 
     920                ssh(:,:,Kbb)=0._wp 
    925921                e3t(:,:,:,Kmm)=e3t_0(:,:,:) 
    926922                e3t(:,:,:,Kbb)=e3t_0(:,:,:) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynatf.F90

    r12529 r12614  
    213213            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    214214               ! Before filtered scale factor at (u/v)-points 
    215                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 
    216                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 
     215               CALL dom_vvl_interpol( ssh(:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 
     216               CALL dom_vvl_interpol( ssh(:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 
    217217               DO_3D_11_11( 1, jpkm1 ) 
    218218                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 
     
    224224               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    225225               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    226                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 
    227                CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 
     226               CALL dom_vvl_interpol( ssh(:,:,Kmm), ze3u_f, 'U' ) 
     227               CALL dom_vvl_interpol( ssh(:,:,Kmm), ze3v_f, 'V' ) 
    228228               DO_3D_11_11( 1, jpkm1 ) 
    229229                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90

    r12529 r12614  
    461461               zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  & 
    462462                  &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     463!!an                  &                        * fmask(ji,jj,jk) 
    463464            END_2D 
    464465         CASE ( np_CME )                           !* Coriolis + metric 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/nemogcm.F90

    r12529 r12614  
    7171   USE stopts         ! Stochastic param.: ??? 
    7272   USE diu_layers     ! diurnal bulk SST and coolskin 
    73    USE step_diu       ! diurnal bulk SST timestepping (called from here if run offline) 
    7473   USE crsini         ! initialise grid coarsening utility 
    7574   USE dia25h         ! 25h mean output 
     
    140139      CALL nemo_init               !==  Initialisations  ==! 
    141140      !                            !-----------------------! 
     141       
    142142#if defined key_agrif 
    143143      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
     
    165165      istp = nit000 
    166166      ! 
    167 #if defined key_c1d 
    168       DO WHILE ( istp <= nitend .AND. nstop == 0 )    !==  C1D time-stepping  ==! 
    169          CALL stp_c1d( istp ) 
     167      !                                               !==  Standard time-stepping  ==! 
     168      ! 
     169      DO WHILE( istp <= nitend .AND. nstop == 0 ) 
     170 
     171         ncom_stp = istp 
     172         IF( ln_timing ) THEN 
     173            zstptiming = MPI_Wtime() 
     174            IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 
     175            IF ( istp ==         nitend ) elapsed_time = zstptiming - elapsed_time 
     176         ENDIF 
     177          
     178         CALL stp        ( istp )  
    170179         istp = istp + 1 
     180 
     181         IF( lwp .AND. ln_timing )   WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 
     182 
    171183      END DO 
    172 #else 
    173       ! 
    174 # if defined key_agrif 
    175       !                                               !==  AGRIF time-stepping  ==! 
    176       CALL Agrif_Regrid() 
    177       ! 
    178       ! Recursive update from highest nested level to lowest: 
    179       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
    180       CALL Agrif_step_child_adj(Agrif_Update_All) 
    181       ! 
    182       DO WHILE( istp <= nitend .AND. nstop == 0 ) 
    183          CALL stp 
    184          istp = istp + 1 
    185       END DO 
    186       ! 
    187       IF( .NOT. Agrif_Root() ) THEN 
    188          CALL Agrif_ParentGrid_To_ChildGrid() 
    189          IF( ln_diaobs )   CALL dia_obs_wri 
    190          IF( ln_timing )   CALL timing_finalize 
    191          CALL Agrif_ChildGrid_To_ParentGrid() 
    192       ENDIF 
    193       ! 
    194 # else 
    195       ! 
    196       IF( .NOT.ln_diurnal_only ) THEN                 !==  Standard time-stepping  ==! 
    197          ! 
    198          DO WHILE( istp <= nitend .AND. nstop == 0 ) 
    199  
    200             ncom_stp = istp 
    201             IF( ln_timing ) THEN 
    202                zstptiming = MPI_Wtime() 
    203                IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 
    204                IF ( istp ==         nitend ) elapsed_time = zstptiming - elapsed_time 
    205             ENDIF 
    206              
    207             CALL stp        ( istp )  
    208             istp = istp + 1 
    209  
    210             IF( lwp .AND. ln_timing )   WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 
    211  
    212          END DO 
    213          ! 
    214       ELSE                                            !==  diurnal SST time-steeping only  ==! 
    215          ! 
    216          DO WHILE( istp <= nitend .AND. nstop == 0 ) 
    217             CALL stp_diurnal( istp )   ! time step only the diurnal SST  
    218             istp = istp + 1 
    219          END DO 
    220          ! 
    221       ENDIF 
    222       ! 
    223 # endif 
    224       ! 
    225 #endif 
    226       ! 
    227       IF( ln_diaobs   )   CALL dia_obs_wri 
    228       ! 
    229       IF( ln_icebergs )   CALL icb_end( nitend ) 
    230  
     184      ! 
     185      ! 
    231186      !                            !------------------------! 
    232187      !                            !==  finalize the run  ==! 
     
    244199      ! 
    245200#if defined key_iomput 
    246                                     CALL xios_finalize  ! end mpp communications with xios 
    247       IF( lk_oasis     )            CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
     201                       CALL xios_finalize  ! end mpp communications with xios 
    248202#else 
    249       IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
    250       ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop      ! end mpp communications 
    251       ENDIF 
     203      IF( lk_mpp   )   CALL mppstop      ! end mpp communications 
    252204#endif 
    253205      ! 
     
    419371      ! 
    420372                           CALL     phy_cst         ! Physical constants 
    421                            CALL     eos_init        ! Equation of state 
    422       IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    423                            CALL     wad_init        ! Wetting and drying options 
     373                            
    424374                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
    425       IF( ln_crs       )   CALL     crs_init(      Nnn )       ! coarsened grid: domain initialization  
     375 
    426376      IF( sn_cfctl%l_prtctl )   & 
    427377         &                 CALL prt_ctl_init        ! Print control 
    428        
    429                            CALL diurnal_sst_bulk_init       ! diurnal sst 
    430       IF( ln_diurnal   )   CALL diurnal_sst_coolskin_init   ! cool skin    
    431       !                             
    432       IF( ln_diurnal_only ) THEN                    ! diurnal only: a subset of the initialisation routines 
    433          CALL  istate_init( Nbb, Nnn, Naa )         ! ocean initial state (Dynamics and tracers) 
    434          CALL     sbc_init( Nbb, Nnn, Naa )         ! Forcings : surface module 
    435          CALL tra_qsr_init                          ! penetrative solar radiation qsr 
    436          IF( ln_diaobs ) THEN                       ! Observation & model comparison 
    437             CALL dia_obs_init( Nnn )                ! Initialize observational data 
    438             CALL dia_obs( nit000 - 1, Nnn )         ! Observation operator for restart 
    439          ENDIF      
    440          IF( lk_asminc )   CALL asm_inc_init( Nbb, Nnn, Nrhs )   ! Assimilation increments 
    441          ! 
    442          RETURN                                       ! end of initialization 
    443       ENDIF 
    444       ! 
    445        
     378 
    446379                           CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
    447380 
    448381      !                                      ! external forcing  
    449382                           CALL    tide_init                     ! tidal harmonics 
     383 
    450384                           CALL     sbc_init( Nbb, Nnn, Naa )    ! surface boundary conditions (including sea-ice) 
    451                            CALL     bdy_init                     ! Open boundaries initialisation 
    452  
    453       !                                      ! Ocean physics 
    454                            CALL zdf_phy_init( Nnn )    ! Vertical physics 
    455                                       
     385                            
     386 
     387      !                                      ! Ocean physics                                     
    456388      !                                         ! Lateral physics 
    457                            CALL ldf_tra_init      ! Lateral ocean tracer physics 
    458                            CALL ldf_eiv_init      ! eddy induced velocity param. 
    459389                           CALL ldf_dyn_init      ! Lateral ocean momentum physics 
    460  
    461       !                                      ! Active tracers 
    462       IF( ln_traqsr    )   CALL tra_qsr_init      ! penetrative solar radiation qsr 
    463                            CALL tra_bbc_init      ! bottom heat flux 
    464                            CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme 
    465                            CALL tra_dmp_init      ! internal tracer damping 
    466                            CALL tra_adv_init      ! horizontal & vertical advection 
    467                            CALL tra_ldf_init      ! lateral mixing 
    468  
     390                            
     391                            
    469392      !                                      ! Dynamics 
    470       IF( lk_c1d       )   CALL dyn_dmp_init         ! internal momentum damping 
    471393                           CALL dyn_adv_init         ! advection (vector or flux form) 
     394 
    472395                           CALL dyn_vor_init         ! vorticity term including Coriolis 
     396 
    473397                           CALL dyn_ldf_init         ! lateral mixing 
    474                            CALL dyn_hpg_init( Nnn )  ! horizontal gradient of Hydrostatic pressure 
     398 
    475399                           CALL dyn_spg_init         ! surface pressure gradient 
    476400 
    477 #if defined key_top 
    478       !                                      ! Passive tracers 
    479                            CALL     trc_init( Nbb, Nnn, Naa ) 
    480 #endif 
    481       IF( l_ldfslp     )   CALL ldf_slp_init    ! slope of lateral mixing 
    482  
    483       !                                      ! Icebergs 
    484                            CALL icb_init( rn_Dt, nit000)   ! initialise icebergs instance 
    485  
    486                                                 ! ice shelf 
    487                            CALL isf_init( Nbb, Nnn, Naa ) 
    488  
    489       !                                      ! Misc. options 
    490                            CALL sto_par_init    ! Stochastic parametrization 
    491       IF( ln_sto_eos   )   CALL sto_pts_init    ! RRandom T/S fluctuations 
    492       
    493401      !                                      ! Diagnostics 
    494402                           CALL     flo_init( Nnn )    ! drifting Floats 
     403                            
    495404      IF( ln_diacfl    )   CALL dia_cfl_init    ! Initialise CFL diagnostics 
    496 !                           CALL dia_ptr_init    ! Poleward TRansports initialization 
    497                            CALL dia_dct_init    ! Sections tranports 
    498                            CALL dia_hsb_init( Nnn )    ! heat content, salt content and volume budgets 
     405 
    499406                           CALL     trd_init( Nnn )    ! Mixed-layer/Vorticity/Integral constraints trends 
    500                            CALL dia_obs_init( Nnn )    ! Initialize observational data 
    501                            CALL dia_25h_init( Nbb )    ! 25h mean  outputs 
    502                            CALL dia_detide_init ! Weights computation for daily detiding of model diagnostics 
    503       IF( ln_diaobs    )   CALL dia_obs( nit000-1, Nnn )   ! Observation operator for restart 
    504                            CALL dia_mlr_init    ! Initialisation of IOM context management for multiple-linear-regression analysis 
    505  
    506       !                                      ! Assimilation increments 
    507       IF( lk_asminc    )   CALL asm_inc_init( Nbb, Nnn, Nrhs )   ! Initialize assimilation increments 
    508       ! 
     407 
     408 
    509409      IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA 
    510410      ! 
    511411      IF( ln_timing    )   CALL timing_stop( 'nemo_init') 
    512412      ! 
     413 
    513414   END SUBROUTINE nemo_init 
    514415 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/phycst.F90

    r12529 r12614  
    1616   !!---------------------------------------------------------------------- 
    1717   USE par_oce          ! ocean parameters 
    18    USE in_out_manager   ! I/O manager 
     18   USE in_out_manager   ! I/O manager  
    1919 
    2020   IMPLICIT NONE 
     
    6666   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos 
    6767   REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi 
     68 
    6869   !!---------------------------------------------------------------------- 
    6970   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    9293      r1_rhos = 1._wp / rhos 
    9394      r1_rcpi = 1._wp / rcpi 
    94  
     95      ! 
     96      rho0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     97      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
     98      ! 
     99      rho0_rcp    = rho0 * rcp  
     100      r1_rho0     = 1._wp / rho0 
     101      r1_rcp      = 1._wp / rcp 
     102      r1_rho0_rcp = 1._wp / rho0_rcp  
     103      ! 
    95104      IF(lwp) THEN 
    96105         WRITE(numout,*) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90

    r12529 r12614  
    244244               ! Horizontal scale factor interpolations 
    245245               ! -------------------------------------- 
    246                CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 
    247                CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 
    248                CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
    249                CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
    250                CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 
     246               CALL dom_vvl_interpol( ssh(:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 
     247               CALL dom_vvl_interpol( ssh(:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 
     248               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     249               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     250               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' ) 
    251251               ! Vertical scale factor interpolations 
    252252               ! ------------------------------------ 
    253                CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
    254                CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
    255                CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
    256                CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
    257                CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
     253               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
     254               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     255               CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     256               CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     257               CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    258258               ! t- and w- points depth 
    259259               ! ---------------------- 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90

    r12529 r12614  
    22   !!====================================================================== 
    33   !!                       ***  MODULE step  *** 
    4    !! Time-stepping   : manager of the ocean, tracer and ice time stepping 
     4   !! Time-stepping   : manager of the shallow water equation time stepping 
    55   !!====================================================================== 
    6    !! History :  OPA  !  1991-03  (G. Madec)  Original code 
    7    !!             -   !  1991-11  (G. Madec) 
    8    !!             -   !  1992-06  (M. Imbard)  add a first output record 
    9    !!             -   !  1996-04  (G. Madec)  introduction of dynspg 
    10    !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer 
    11    !!            8.0  !  1997-06  (G. Madec)  new architecture of call 
    12    !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface 
    13    !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit 
    14    !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions 
    15    !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking 
    16    !!             -   !  2004-08  (C. Talandier) New trends organization 
    17    !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme 
    18    !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls 
    19    !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation 
    20    !!             -   !  2006-07  (S. Masson)  restart using iom 
    21    !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate 
    22    !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl 
    23    !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    24    !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 
    25    !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal 
    26    !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs 
    27    !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state 
    28    !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 
    29    !!            3.7  !  2014-10  (G. Madec)  LDF simplication  
    30    !!             -   !  2014-12  (G. Madec) remove KPP scheme 
    31    !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface) 
    32    !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy) 
    33    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 
    34    !!---------------------------------------------------------------------- 
    35  
    36    !!---------------------------------------------------------------------- 
    37    !!   stp             : OPA system time-stepping 
     6   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2 
     7   !!---------------------------------------------------------------------- 
     8 
     9   !!---------------------------------------------------------------------- 
     10   !!   stp             : Shallow Water time-stepping 
    3811   !!---------------------------------------------------------------------- 
    3912   USE step_oce         ! time stepping definition modules 
     13   USE phycst           ! physical constants 
     14   USE usrdef_nam 
    4015   ! 
    41    USE iom              ! xIOs server 
     16   USE iom              ! xIOs server  
    4217 
    4318   IMPLICIT NONE 
     
    4520 
    4621   PUBLIC   stp   ! called by nemogcm.F90 
    47  
     22    
    4823   !!---------------------------------------------------------------------- 
    4924   !! time level indices 
    5025   !!---------------------------------------------------------------------- 
    51    INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init 
    52  
     26   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init 
     27       
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    5330   !!---------------------------------------------------------------------- 
    5431   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6845      !!                     ***  ROUTINE stp  *** 
    6946      !! 
    70       !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.) 
    71       !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.) 
    72       !!              - Time stepping of TRC  (passive tracer eqs.) 
     47      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.) 
    7348      !! 
    74       !! ** Method  : -1- Update forcings and data 
    75       !!              -2- Update ocean physics 
    76       !!              -3- Compute the t and s trends 
    77       !!              -4- Update t and s 
    78       !!              -5- Compute the momentum trends 
    79       !!              -6- Update the horizontal velocity 
    80       !!              -7- Compute the diagnostics variables (rd,N2, hdiv,w) 
    81       !!              -8- Outputs and diagnostics 
     49      !! ** Method  : -1- Update forcings 
     50      !!              -2- Update the ssh at Naa 
     51      !!              -3- Compute the momentum trends (Nrhs) 
     52      !!              -4- Update the horizontal velocity 
     53      !!              -5- Apply Asselin time filter to uu,vv,ssh 
     54      !!              -6- Outputs and diagnostics 
    8255      !!---------------------------------------------------------------------- 
    8356      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     
    8558!!gm kcall can be removed, I guess 
    8659      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt) 
     60      REAL(wp)::   z1_2rho0     ! local scalars 
     61       
     62      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars 
     63      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      - 
     64      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 
    8765      !! --------------------------------------------------------------------- 
    8866#if defined key_agrif 
     
    10583      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    10684      ! 
    107       IF( l_1st_euler ) THEN   
    108          ! start or restart with Euler 1st time-step 
    109          rDt =  rn_Dt    
     85      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step 
     86         rDt   =  rn_Dt    
    11087         r1_Dt = 1._wp / rDt 
    11188      ENDIF 
     89 
     90      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero 
     91 
    11292      ! 
    11393      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    132112      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    133113      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries 
    134       IF( ln_isf     )   CALL isf_stp ( kstp, Nnn ) 
    135                          CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice) 
    136  
    137       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    138       ! Update stochastic parameters and random T/S fluctuations 
    139       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    140       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters 
    141       IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations 
     114                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition 
    142115 
    143116      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    144117      ! Ocean physics update 
    145118      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    146       !  THERMODYNAMICS 
    147                          CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
    148                          CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
    149                          CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
    150                          CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
    151  
    152       !  VERTICAL PHYSICS 
    153                          CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
    154  
    155119      !  LATERAL  PHYSICS 
    156       ! 
    157       IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
    158                          CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
    159  
    160          IF( ln_zps .AND. .NOT. ln_isfcav)                                    & 
    161             &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    162             &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level 
    163  
    164          IF( ln_zps .AND.       ln_isfcav)                                                & 
    165             &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    166             &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level 
    167          IF( ln_traldf_triad ) THEN  
    168                          CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator 
    169          ELSE      
    170                          CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator 
    171          ENDIF 
    172       ENDIF 
    173120      !                                                                        ! eddy diffusivity coeff. 
    174       IF( l_ldftra_time .OR. l_ldfeiv_time )   CALL ldf_tra( kstp, Nbb, Nnn )  !       and/or eiv coeff. 
    175121      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff.  
    176122 
     
    180126 
    181127                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor) 
     128 
    182129      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors  
    183                             CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity  
    184       IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning 
    185                             CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation 
    186                              
    187                              
     130                                                        
    188131                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    189132                         vv(:,:,:,Nrhs) = 0._wp 
    190133 
    191       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    192                &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment 
    193134      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
     135 
    194136#if defined key_agrif 
    195137      IF(.NOT. Agrif_Root())  &  
     
    197139#endif 
    198140                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
     141  
    199142                         CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
     143  
    200144                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    201       IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
    202                          CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure 
    203                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    204  
    205                                                       ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 
    206       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
    207                             CALL div_hor       ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case) 
    208          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component) 
    209       ENDIF 
    210                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
    211       IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated 
    212                             CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )             ! now cross-level velocity  
    213          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning 
    214       ENDIF 
    215        
    216  
    217       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    218       ! cool skin 
    219       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    220       IF ( ln_diurnal )  CALL diurnal_layers( kstp ) 
    221        
    222       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    223       ! diagnostics and outputs 
    224       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    225       IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats 
    226       IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics 
    227                          CALL dia_hth   ( kstp,      Nnn )      ! Thermocline depth (20 degres isotherm depth) 
    228       IF( ln_diadct  )   CALL dia_dct   ( kstp,      Nnn )      ! Transports 
    229                          CALL dia_ar5   ( kstp,      Nnn )      ! ar5 diag 
    230                          CALL dia_ptr   ( kstp,      Nnn )      ! Poleward adv/ldf TRansports diagnostics 
    231                          CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs 
    232       IF( ln_crs     )   CALL crs_fld   ( kstp,      Nnn )      ! ocean model: online field coarsening & output 
    233       IF( lk_diadetide ) CALL dia_detide( kstp )                ! Weights computation for daily detiding of model diagnostics 
    234       IF( lk_diamlr  )   CALL dia_mlr                           ! Update time used in multiple-linear-regression analysis 
    235        
    236 #if defined key_top 
    237       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    238       ! Passive Tracer Model 
    239       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    240                          CALL trc_stp       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! time-stepping 
    241 #endif 
    242  
    243       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    244       ! Active tracers                               
    245       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    246                          ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero 
    247  
    248       IF(  lk_asminc .AND. ln_asmiau .AND. & 
    249          & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment 
    250                          CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition 
    251       IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr 
    252       IF( ln_isf     )   CALL tra_isf    ( kstp,      Nnn, ts, Nrhs )  ! ice shelf heat flux 
    253       IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux 
    254       IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme 
    255       IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends 
    256       IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends 
    257 #if defined key_agrif 
    258       IF(.NOT. Agrif_Root())  &  
    259                &         CALL Agrif_Sponge_tra        ! tracers sponge 
    260 #endif 
    261                          CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection ==> RHS 
    262       IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS 
    263       IF( lrst_oce .AND. ln_zdfosm ) & 
    264            &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts 
    265                          CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing 
    266  
    267                          CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vertical mixing and after tracer fields 
    268       IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection 
     145 
     146!!an - calcul du gradient de pression horizontal (explicit) 
     147      DO_3D_00_00( 1, jpkm1 ) 
     148         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
     149         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
     150      END_3D 
     151      ! 
     152      ! add wind stress forcing and layer linear friction to the RHS  
     153      z1_2rho0 = 0.5_wp * r1_rho0 
     154      DO_3D_00_00(1,jpkm1) 
     155         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   & 
     156            &                                  - rn_rfr * uu(ji,jj,jk,Nbb) 
     157         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   & 
     158            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)   
     159      END_3D    
     160!!an          
     161   
     162      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     163      ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3  
     164      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     165       
     166!!an futur module dyn_nxt (a la place de dyn_atf) 
     167       
     168      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity 
     169         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter) 
     170            DO_3D_00_00(1,jpkm1) 
     171               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     172               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     173             END_3D           
     174         ELSE                             ! Leap Frog time stepping + Asselin filter          
     175            DO_3D_11_11(1,jpkm1) 
     176               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     177               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     178               !                                                                  ! Asselin time filter on u,v (Nnn) 
     179               uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua) 
     180               vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva) 
     181               ! 
     182               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) ) 
     183               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) ) 
     184               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) ) 
     185               ! 
     186               e3u(ji,jj,jk,Nnn) = ze3u_tf 
     187               e3v(ji,jj,jk,Nnn) = ze3v_tf 
     188               e3t(ji,jj,jk,Nnn) = ze3t_tf 
     189               !               
     190               uu(ji,jj,jk,Naa) = zua 
     191               vv(ji,jj,jk,Naa) = zva 
     192             END_3D    
     193         ENDIF 
     194         ! 
     195      ELSE                          ! flux form : applied on thickness weighted velocity 
     196         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter) 
     197            DO_3D_00_00(1,jpkm1) 
     198               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
     199               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
     200               !                                                ! LF time stepping 
     201               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     202               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     203               !                                                 
     204               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     205               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)      
     206             END_3D           
     207         ELSE                             ! Leap Frog time stepping + Asselin filter          
     208            DO_3D_11_11(1,jpkm1) 
     209               zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) 
     210               zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn) 
     211               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
     212               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
     213               !                                                ! LF time stepping 
     214               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     215               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     216               !                                                ! Asselin time filter on e3u/v/t 
     217               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) )  
     218               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) ) 
     219               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) )  
     220               !                                                ! Asselin time filter on u,v (Nnn) 
     221               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_tf 
     222               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_tf            
     223               ! 
     224               e3u(ji,jj,jk,Nnn) = ze3u_tf 
     225               e3v(ji,jj,jk,Nnn) = ze3v_tf 
     226               e3t(ji,jj,jk,Nnn) = ze3t_tf 
     227               ! 
     228               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     229               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)      
     230             END_3D    
     231         ENDIF 
     232      ENDIF 
     233       
     234      CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries 
     235         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )      
     236 
     237!!an          
    269238 
    270239      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    271240      ! Set boundary conditions, time filter and swap time levels 
    272241      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    273 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap  
    274 !!    (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.  
    275 !!    If so:  
    276 !!    (i) no need to call agrif update at initialization time 
    277 !!    (ii) no need to update "before" fields  
    278 !! 
    279 !!    Apart from creating new tra_swp/dyn_swp routines, this however:  
    280 !!    (i) makes boundary conditions at initialization time computed from updated fields which is not the case between  
    281 !!    two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",  
    282 !!    e.g. a shift of the feedback interface inside child domain.  
    283 !!    (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 
    284 !!    place. 
    285 !!  
    286 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 
    287                          CALL tra_atf       ( kstp, Nbb, Nnn, Naa, ts )                      ! time filtering of "now" tracer arrays 
    288                          CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
    289                          CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height 
    290       ! 
     242 
     243!!an TO BE ADDED : dyn_nxt 
     244!!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors 
     245!!an TO BE ADDED : a simplifier 
     246!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height 
     247  
     248      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps 
     249         !                                                  ! filtering "now" field 
     250         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 
     251      ENDIF 
     252       
     253!!an  
     254 
     255 
    291256      ! Swap time levels 
    292257      Nrhs = Nbb 
     
    295260      Naa = Nrhs 
    296261      ! 
    297       IF(.NOT.ln_linssh) CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors 
    298       ! 
    299       IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics 
    300  
    301 !!gm : This does not only concern the dynamics ==>>> add a new title 
    302 !!gm2: why ouput restart before AGRIF update? 
    303 !! 
    304 !!jc: That would be better, but see comment above 
    305 !! 
     262                         CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors 
     263 
     264      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     265      ! diagnostics and outputs 
     266      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     267      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats 
     268      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics 
     269     
     270                        CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs 
     271 
     272      ! 
    306273      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file 
    307       IF( ln_sto_eos )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
     274           
    308275 
    309276#if defined key_agrif 
     
    324291      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    325292                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic ) 
     293                
    326294                          
    327295      IF( kstp == nit000 ) THEN                          ! 1st time step only 
     
    331299      ENDIF 
    332300 
    333       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    334       ! Coupled mode 
    335       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    336 !!gm why lk_oasis and not lk_cpl ???? 
    337       IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges 
    338301      ! 
    339302#if defined key_iomput 
     
    341304                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    342305                      IF(lrxios) CALL iom_context_finalize(      crxios_context          ) 
    343          IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !  
    344306      ENDIF 
    345307#endif 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90

    r12529 r12614  
    6767      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
    6868      REAL(wp)               ::   zzz                 ! local real  
    69       REAL(wp), DIMENSION(9) ::   zmax 
     69      REAL(wp), DIMENSION(3) ::   zmax 
    7070      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
    7171      CHARACTER(len=20) :: clname 
     
    9191            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 
    9292            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu   ) 
    93             istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  ) 
    94             istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  ) 
    95             istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  ) 
    96             istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  ) 
    97             IF( ln_zad_Aimp ) THEN 
    98                istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  ) 
    99                istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
    100             ENDIF 
    10193            istatus = NF90_ENDDEF(idrun) 
    102             zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use 
    10394         ENDIF 
    10495      ENDIF 
     
    114105         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
    115106      ELSE 
    116          zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max 
     107         zmax(1) = MINVAL( e3t(:,:,1,Kmm)  )                                         ! ssh min 
    117108      ENDIF 
    118109      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
    119       zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
    120       zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
    121       zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
    122       zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
    123       zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
    124       IF( ln_zad_Aimp ) THEN 
    125          zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
    126          zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
    127       ENDIF 
     110      zmax(3) = REAL( nstop , wp )                                            ! stop indicator 
    128111      ! 
    129112      IF( ll_colruns ) THEN 
    130113         CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    131          nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains 
     114         nstop = NINT( zmax(3) )                 ! nstop indicator sheared among all local domains 
    132115      ENDIF 
    133116      !                                   !==  run statistics  ==!   ("run.stat" files) 
    134117      IF( ll_wrtruns ) THEN 
    135          WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
     118         WRITE(numrun,9500) kt, zmax(1), zmax(2) 
    136119         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 
    137120         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) ) 
    138          istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) ) 
    139          istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) ) 
    140          istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) ) 
    141          istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) ) 
    142          IF( ln_zad_Aimp ) THEN 
    143             istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) ) 
    144             istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) ) 
    145          ENDIF 
    146121         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
    147122         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
     
    149124      !                                   !==  error handling  ==! 
    150125      IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges 
    151          &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
     126         &  zmax(1) <    0._wp .OR.   &                    ! negative sea surface height  
    152127         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    153          &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
    154          &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
    155          &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
    156          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
     128         &  ISNAN( zmax(1) + zmax(2) ) ) ) THEN            ! NaN encounter in the tests 
    157129         IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 
    158130            ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 
    159131            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
    160132            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
    161             CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 
    162             CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 
    163133         ELSE 
    164134            ! find local min and max locations 
    165135            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
    166136            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    167             is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    168             is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    169137         ENDIF 
    170138          
    171          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
     139         WRITE(ctmp1,*) ' stp_ctl: (e3t0) ssh < 0 m  or  |U| > 10 m/s  or  NaN encounter in the tests' 
    172140         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
    173141         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
    174          WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 
    175          WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3) 
    176          WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file' 
     142         WRITE(ctmp4,*) '      ===> output of last computed fields in output.abort.nc file' 
    177143          
    178144         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
     
    180146         IF( .NOT. sn_cfctl%l_glochk ) THEN 
    181147            WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 
    182             CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 ) 
     148            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4 ) 
    183149         ELSE 
    184             CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' ) 
     150            CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4 ) 
    185151         ENDIF 
    186152 
     
    189155      ENDIF 
    190156      ! 
    191 9100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5) 
     1579100  FORMAT (' kt=',i8,'   |ssh| min: ',1pg11.4,', at  i j  : ',2i5) 
    1921589200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
    193 9300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5) 
    194 9400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5) 
    195 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
     1599500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16) 
    196160      ! 
    197161   END SUBROUTINE stp_ctl 
Note: See TracChangeset for help on using the changeset viewer.