Changeset 13416


Ignore:
Timestamp:
2020-08-20T12:10:55+02:00 (5 months ago)
Author:
gm
Message:

First dev of BVP, cond 3.4, gridF fixes

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater
Files:
2 added
11 edited
5 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/field_def_nemo-oce.xml

    r13047 r13416  
    645645      <!-- T grid --> 
    646646      <field_group id="grid_T" grid_ref="grid_T_2D">   
    647         <field id="pry"         long_name="porosity parameter"      standard_name="porosity_phi"         unit="#"    /> 
    648         <field id="pmb"         long_name="permeability parameter"  standard_name="permeability_sigma"   unit="1/s" /> 
     647        <field id="rpo"         long_name="porosity parameter"      standard_name="porosity_phi"         unit=""    /> 
     648        <field id="bmpt"        long_name="permeability parameter T"  standard_name="permeability_sigma_t"   unit="1/s" /> 
     649      </field_group> 
     650       
     651      <!-- U grid --> 
     652      <field_group id="grid_U" grid_ref="grid_U_2D">   
     653        <field id="bmpu"         long_name="permeability parameter U facets"  standard_name="permeability_sigma_u"   unit="1/s" /> 
     654      </field_group> 
     655 
     656      <!-- V grid -->       
     657      <field_group id="grid_V" grid_ref="grid_V_2D">   
     658        <field id="bmpv"         long_name="permeability parameter V facets"  standard_name="permeability_sigma_v"   unit="1/s" /> 
    649659      </field_group> 
    650660       
     
    659669        <field id="Ens"          long_name="enstrophy"                         standard_name="enstrophy"              unit="1/m2.s2"               /> 
    660670      </field_group> 
    661        
    662        
     671  
     672      <field_group id="grid_F" grid_ref="grid_F_3D">   
     673        <field id="e3f_0"        long_name="F-cell thickness"                  standard_name="cell_thickness"        unit="m"  /> 
     674        <field id="e3f"          long_name="F-cell thickness"                  standard_name="cell_thickness"        unit="m"  /> 
     675      </field_group>      
     676  
    663677      <!-- AGRIF sponge --> 
    664678      <field id="agrif_spf"    long_name=" AGRIF f-sponge coefficient"   unit=" " /> 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/file_def_nemo-oce.xml

    r12667 r13416  
    2222       
    2323      <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
    24           <field field_ref="ht"           name="ht"       operation = "instant" /> 
    25           <field field_ref="sKE"          name="sKE"      operation = "instant" /> 
    26           <field field_ref="ssh"          name="sossheig" operation = "instant" /> 
    27           <field field_ref="wetdep"       name="hswe_wd"  operation = "instant" /> 
     24          <field field_ref="e3t"          name="e3t"       operation = "instant" /> 
     25          <field field_ref="e3t_0"        name="e3t_0"     operation = "instant" /> 
     26          <field field_ref="ht"           name="ht"        operation = "instant" /> 
     27          <field field_ref="sKE"          name="sKE"       operation = "instant" /> 
     28          <field field_ref="ssh"          name="sossheig"  operation = "instant" /> 
     29          <field field_ref="wetdep"       name="hswe_wd"   operation = "instant" /> 
     30          <field field_ref="rpo"          name="rpo"       operation = "instant" /> 
     31          <field field_ref="bmpt"         name="bmpt"      operation = "instant" /> 
    2832        </file> 
    2933 
    3034        <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" > 
    31           <field field_ref="hu"           name="hu"       operation = "instant" /> 
    32           <field field_ref="ssu"          name="ssu"      operation = "instant"  /> 
    33           <field field_ref="utau"         name="sozotaux" operation = "instant"  /> 
     35          <field field_ref="e3u"          name="e3u"        operation = "instant" /> 
     36          <field field_ref="e3u_0"        name="e3u_0"      operation = "instant" /> 
     37          <field field_ref="hu"           name="hu"         operation = "instant" /> 
     38          <field field_ref="ssu"          name="ssu"        operation = "instant"  /> 
     39          <field field_ref="utau"         name="sozotaux"   operation = "instant"  /> 
     40          <field field_ref="bmpt"         name="bmpv"       operation = "instant" /> 
    3441        </file> 
    3542 
    3643        <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" > 
    37           <field field_ref="hv"           name="hv"       operation = "instant" /> 
    38           <field field_ref="ssv"          name="ssv"      operation = "instant"  /> 
    39           <field field_ref="vtau"         name="sometauy" operation = "instant"  /> 
     44          <field field_ref="e3v"          name="e3v"        operation = "instant" /> 
     45          <field field_ref="e3v_0"        name="e3v_0"      operation = "instant" /> 
     46          <field field_ref="hv"           name="hv"         operation = "instant" /> 
     47          <field field_ref="ssv"          name="ssv"        operation = "instant"  /> 
     48          <field field_ref="vtau"         name="sometauy"   operation = "instant"  /> 
     49          <field field_ref="bmpv"         name="bmpv"       operation = "instant" /> 
    4050        </file> 
    4151         
    4252        <file id="file5" name_suffix="_grid_F" description="ocean F grid variables put on T grid" > 
     53          <field field_ref="e3f"            name="e3f"       operation = "instant" /> 
     54          <field field_ref="e3f_0"          name="e3f_0"     operation = "instant" /> 
    4355          <field field_ref="hf"             name="hf"        operation = "instant" /> 
    4456          <field field_ref="sKEf"           name="sKEf"      operation = "instant" /> 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg

    r13005 r13416  
    6262   !                            !       = 3   symmetric laplacian (cartesian) 
    6363   ln_dynldf_lap_PM  =.false.   ! if True - apply the P.Marchand boundary condition on the laplacian viscosity 
     64   ! 
     65   rn_abp = 1.e-2     ! alpha boundary parameter                                       [-] 
     66   nn_cnp = 1         ! number of cell on which is smoothed the porosity (phi)         [-] 
     67   rn_fsp = 0.        ! friction parameter 1/epsilon of the permeability                 [1/s] 
    6468   ! 
    6569/ 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/diawri.F90

    r13151 r13416  
    7575   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file 
    7676   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
     77   INTEGER ::   nid_F, nz_F, nh_F, ndim_F, ndim_hF   ! grid_F file 
    7778   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
    7879   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7980   INTEGER ::   ndex(1)                              ! ??? 
    80    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV, ndex_hF 
    8182   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    82    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
     83   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V, ndex_F 
    8384   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    8485 
     
    138139      CALL iom_put("e3u_0", e3u_0(:,:,:) ) 
    139140      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
     141      CALL iom_put("e3f_0", e3f_0(:,:,:) ) 
    140142      ! 
    141143!!st13 
     
    173175         CALL iom_put( "ssh" , ssh(:,:,Kmm) )                          ! sea surface height 
    174176      ENDIF 
    175  
    176 !!an  
    177         IF( iom_use("ht") )   &                  ! water column at t-point 
     177      ! 
     178# if defined key_bvp 
     179      IF( iom_use("rpo") )    CALL iom_put( "rpo", rpo(:,:,1) ) 
     180      IF( iom_use("bmpt") )   CALL iom_put( "bmpt", bmpt(:,:,1) ) 
     181      IF( iom_use("bmpu") )   CALL iom_put( "bmpu", bmpu(:,:,1) ) 
     182      IF( iom_use("bmpv") )   CALL iom_put( "bmpv", bmpv(:,:,1) ) 
     183#endif  
     184      ! 
     185      IF( iom_use("ht") )   &                  ! water column at t-point 
    178186         CALL iom_put( "ht" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    179187      ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/domain.F90

    r13151 r13416  
    154154      CALL dom_zgr( ik_top, ik_bot )    ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices) 
    155155 
    156       CALL dom_msk( ik_top, ik_bot )    ! Masks 
     156      CALL dom_msk( ik_top, ik_bot )    ! Masks  
     157                                        ! Penalisation fields 
     158      ! 
     159# if defined key_bvp 
     160      !  Multiply the scale factors by the porosity field 
     161      ! ------------------------------------------------- 
     162      e3t_0(:,:,:) = rpo(:,:,:) * e3t_0(:,:,:) 
     163      ! phiUe3U = mean_i(phiTe3T) 
     164      DO_3D_00_00(1,jpkm1) 
     165         e3u_0(ji,jj,jk) = 0.5_wp  * ( e3t_0(ji,jj,jk) + e3t_0(ji+1,jj,jk) )  
     166         e3v_0(ji,jj,jk) = 0.5_wp  * ( e3t_0(ji,jj,jk) + e3t_0(ji,jj+1,jk) )   
     167         e3f_0(ji,jj,jk) = 0.25_wp * ( e3t_0(ji,jj+1,jk) + e3t_0(ji+1,jj+1,jk)   & 
     168            &                        + e3t_0(ji,jj  ,jk) + e3t_0(ji+1,jj  ,jk)   ) 
     169      END_3D 
     170      ! copy the communication layer (kinf od periodic over the basin) 
     171      CALL lbc_lnk_multi( 'domain', e3t_0, 'T', 1., e3u_0, 'U', 1.,                      & 
     172              &                     e3v_0, 'V', 1., e3f_0, 'F', 1., kfillmode=jpfillcopy ) 
     173      !     
     174      ! why not (ca se fait une fois, ca ne coute pas si cher) 
     175      ! ou, WHERE( tmask(:,:,:) /= 0._wp )   e3t_0(:,:,:) = pry(:,:,:) * e3t_0(:,:,:) 
     176      ! phiU = mean_i(phiT) 
     177      !DO_3D_00_00 (1,jpkm1)                            
     178      !   e3u_0(ji,jj,jk) = 0.5_wp * ( pry(ji,jj,jk) + pry(ji+1,jj,jk) )  * e3u_0(ji,jj,jk) 
     179      !   e3v_0(ji,jj,jk) = 0.5_wp * ( pry(ji,jj,jk) + pry(ji,jj+1,jk) )  * e3v_0(ji,jj,jk) 
     180      !   e3f_0(ji,jj,jk) = 0.25_wp * ( pry(ji,jj+1,jk) + pry(ji+1,jj+1,jk)   & 
     181      !      &                          pry(ji,jj  ,jk) + pry(ji+1,jj  ,jk)   ) * e3f_0(ji,jj,jk) 
     182      !END_3D 
     183      !CALL lbc_lnk_multi( 'domain', e3u_0, 'U', 1., e3v_0, 'V', 1., e3f_0, 'F', 1., kfillmode=jfillcopy )     
     184#endif  
    157185      ! 
    158186      ht_0(:,:) = 0._wp  ! Reference ocean thickness 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/dommsk.F90

    r13151 r13416  
    3232   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3333   USE lib_mpp        ! Massively Parallel Processing library 
     34   ! 
     35   USE usrdef_nam , ONLY : rn_abp, r1_abp, rn_fsp 
    3436 
    3537   IMPLICIT NONE 
     
    9294      INTEGER  ::   iktop, ikbot   !   -       - 
    9395      INTEGER  ::   ios, inum 
    94       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwf   ! 2D workspace 
     96      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zwf   ! 2D workspace 
     97      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3   ! 3D workspace 
    9598      !! 
    9699      NAMELIST/namlbc/ rn_shlat, ln_vorlat 
     
    139142      END_2D 
    140143      ! 
     144# if defined key_bvp 
     145      !  Add a land layer around the domain of porosity alpha  
     146      ! ----------------------------------------------------- 
     147!!an to be added : lissage sur nn_cnp dx  
     148      WRITE(numout,*) 'dommsk rpo entering' 
     149      !zw3   (:,:,:) = tmask(:,:,:) 
     150      zw3 = tmask 
     151      rpo   (:,:,:) = 1._wp  
     152      r1_rpo(:,:,:) = 1._wp 
     153      bmpt  (:,:,:) = 0._wp 
     154      ! alternative pour hf dans les angles 
     155      WHERE(tmask == 0._wp) 
     156         rpo   (:,:,:) = rn_abp             ! used on T points 
     157         r1_rpo(:,:,:) = r1_abp             !  
     158      END WHERE 
     159      WRITE(numout,*) 'dommsk rpo loop' 
     160      ! futur : boucle 3D et test 2D sur les tmask(ji,jj,jk) pour choper la bathy 
     161      DO_2D_00_00                            
     162         IF( tmask(ji,jj,1) < tmask(ji+1,jj+1,1) .OR.   &    
     163          &  tmask(ji,jj,1) < tmask(ji+1,jj-1,1) .OR.   & 
     164          &  tmask(ji,jj,1) < tmask(ji-1,jj-1,1) .OR.   & 
     165          &  tmask(ji,jj,1) < tmask(ji-1,jj+1,1)        )   THEN   
     166          ! plus simple : tmask(ji) + tmask(ji+1)= 0.5  idem pour jj en .OR. 
     167          ! mais j'ai pas mon coin 
     168          ! pour smooth, boucle un peu similaire, si pénalisé alors son copain l'est aussi mais un peu moins 
     169          ! peut etre sur pry 
     170          ! to be added : pry = 0 dans la terre et r1_pry vaut 1 pour ne pas avoir de 1/0 dans certain cas 
     171          ! non, car e3t0 vaudra 0 alors 
     172            !                                           ! flat bottom  
     173            zw3   (ji,jj, 1:jpkm1) = 1._wp 
     174            !rpo   (ji,jj, 1:jpkm1) = rn_abp             ! used on T points 
     175            !r1_rpo(ji,jj, 1:jpkm1) = r1_abp             !  
     176            bmpt  (ji,jj, 1:jpkm1) = rn_fsp             !  
     177         ENDIF 
     178      END_2D 
     179      tmask = zw3 
     180      CALL lbc_lnk_multi( 'dommsk', rpo,  'T', 1._wp, r1_rpo  , 'T', 1._wp,                     & 
     181              &                     bmpt, 'T', 1._wp,                       kfillmode=jpfillcopy )                  
     182      WRITE(numout,*) 'dommsk rpo done' 
     183#endif 
     184      ! 
    141185      ! the following call is mandatory 
    142186      ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...)   
     
    157201         END_3D 
    158202      ENDIF 
    159           
     203      ! 
    160204      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask) 
    161205      ! ---------------------------------------- 
     
    174218      END DO 
    175219      CALL lbc_lnk_multi( 'dommsk', umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions 
    176   
     220 
     221# if defined key_bvp 
     222      !  Permeability applied on the inner boundary layer  
     223      ! -------------------------------------------- 
     224      WRITE(numout,*) 'dommsk bmp loop' 
     225      bmpu(:,:,:) = 0._wp 
     226      bmpv(:,:,:) = 0._wp 
     227      DO_2D_00_00                            
     228         IF(    ( umask(ji-1,jj+1,1) + umask(ji+1,jj+1,1)   &        ! U points 
     229            &   + umask(ji-1,jj-1,1) + umask(ji+1,jj-1,1)   ) == 3._wp) THEN    
     230            ! 
     231            bmpu(ji,jj, 1:jpkm1) = rn_fsp 
     232            ! 
     233         ENDIF 
     234         IF(    ( vmask(ji-1,jj+1,1) + vmask(ji+1,jj+1,1)   &        ! V points 
     235            &   + vmask(ji-1,jj-1,1) + vmask(ji+1,jj-1,1)   ) == 3._wp) THEN    
     236            ! 
     237            bmpv(ji,jj, 1:jpkm1) = rn_fsp 
     238            ! 
     239         ENDIF 
     240      END_2D 
     241      WRITE(numout,*) 'dommsk bmp done' 
     242      CALL lbc_lnk_multi( 'dommsk', bmpu,  'U', 1._wp, bmpv  , 'V', 1._wp, kfillmode=jpfillcopy )    
     243#endif 
     244        
    177245      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask) 
    178246      !----------------------------------------- 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/dynldf_lap_blp.F90

    r13151 r13416  
    2020   USE in_out_manager ! I/O manager 
    2121   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    22  
     22   USE lib_mpp        ! MPP library 
     23   ! 
     24   USE usrdef_nam , ONLY : nn_dynldf_lap_typ    ! use laplacian parameter 
     25   ! 
    2326   IMPLICIT NONE 
    2427   PRIVATE 
     
    3134   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian) 
    3235    
    33    INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist) 
     36   !INTEGER, PUBLIC, PARAMETER ::   nn_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist) 
    3437!!anSYM 
    3538   !! * Substitutions 
     
    7073      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
    7174      REAL(wp), DIMENSION(jpi,jpj) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms 
     75      ! 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t,ze3u,ze3v, ze3f     ! not penalised scale factors 
    7277      !!---------------------------------------------------------------------- 
    7378      ! 
     
    7984         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 
    8085         WRITE(numout,*) '~~~~~~~ ' 
    81          WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ 
    82          SELECT CASE( ln_dynldf_lap_typ )             ! print the choice of operator 
     86         WRITE(numout,*) '                  nn_dynldf_lap_typ = ', nn_dynldf_lap_typ 
     87         SELECT CASE( nn_dynldf_lap_typ )             ! print the choice of operator 
    8388         CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian' 
    8489         CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)' 
    85          CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)' 
     90         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (cartesian form)' 
    8691         END SELECT 
    8792      ENDIF 
     
    9196      ENDIF 
    9297      ! 
    93       SELECT CASE( ln_dynldf_lap_typ )   
     98!!an  Kbb because grad et LeapFrog explicit 
     99!!      ze3t(:,:,:) = e3t(:,:,:,Kbb) 
     100!!      ze3u(:,:,:) = e3u(:,:,:,Kmm) 
     101!!      ze3v(:,:,:) = e3v(:,:,:,Kmm) 
     102!!      ze3f(:,:,:) = e3f(:,:,:) 
     103!!# if defined key_bvp 
     104!!      !   gradient and divergence not penalised 
     105!!      ze3t(:,:,:) = ze3t(:,:,:) * r1_rpo(:,:,:)   
     106!!      ze3u(:,:,:) = ze3u(:,:,:) * r1_rpo(:,:,:)   
     107!!      ze3v(:,:,:) = ze3v(:,:,:) * r1_rpo(:,:,:)   
     108!!      !ze3f(:,:,:) = ze3f(:,:,:) * r1_rpo(:,:,:)    ! unused 
     109!!#endif 
     110      SELECT CASE( nn_dynldf_lap_typ )   
    94111         !               
    95112         CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==! 
     
    99116               DO_2D_01_01 
    100117               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    101 !!gm open question here : e3f  at before or now ?    probably now...  
    102118!!gm note that ahmf has already been multiplied by fmask 
    103119            zcur(ji-1,jj-1) =  & 
     
    193209            !   
    194210         CASE DEFAULT                                     ! error 
    195             CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ'  ) 
     211            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ'  ) 
    196212         END SELECT 
    197213         ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/sshwzv.F90

    r13151 r13416  
    7777      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height 
    7878      ! 
    79       INTEGER  ::   jk      ! dummy loop index 
     79      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    8080      REAL(wp) ::   zcoef   ! local scalar 
    8181      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t    ! penalised scale factors 
    8283      !!---------------------------------------------------------------------- 
    8384      ! 
     
    9899         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv ) 
    99100      ENDIF 
    100  
     101      ! 
     102      DO_3D_00_00( 1, jpkm1 ) 
     103         ze3t(ji,jj,jk) = e3t(ji,jj,jk,Kmm) 
     104      END_3D 
     105# if defined key_bvp 
     106      !   divergence not penalised 
     107      ze3t(:,:,:) = ze3t(:,:,:) * r1_rpo(:,:,:)   
     108#endif 
    101109      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence 
    102110      ! 
    103111      zhdiv(:,:) = 0._wp 
    104112      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    105         zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
     113        zhdiv(:,:) = zhdiv(:,:) + ze3t(:,:,jk) * hdiv(:,:,jk) 
    106114      END DO 
    107115      !                                                ! Sea surface elevation time stepping 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90

    r13005 r13416  
    4242   !                                                  ! = 3   symmetric laplacian (cartesian) 
    4343   LOGICAL , PUBLIC ::   ln_dynldf_lap_PM             ! if True - apply the P.Marchand boundary condition on the laplacian 
     44   !                              !!* penalisation *!! 
     45   REAL(wp), PUBLIC ::   rn_abp             ! alpha boundary parameter                                       [-] 
     46   INTEGER , PUBLIC ::   nn_cnp             ! number of cell on which is smoothed the porosity (phi)         [-] 
     47   REAL(wp), PUBLIC ::   rn_fsp             ! friction parameter 1/epsilon of the permeability               [1/s] 
     48   ! 
     49   REAL(wp), PUBLIC ::   r1_abp             ! inverse alpha boundary parameter                            [-] 
    4450   ! 
    4551   !!---------------------------------------------------------------------- 
     
    7480         &                 rn_f0 ,rn_beta,                      &   ! coriolis parameter 
    7581         &                 rn_modified_grav, rn_rfr , rn_tau,   &   ! reduced gravity, friction, wind 
    76          &                 nn_dynldf_lap_typ, ln_dynldf_lap_PM      ! temporary parameter 
     82         &                 nn_dynldf_lap_typ, ln_dynldf_lap_PM, &   ! temporary parameter 
     83         &                 rn_abp, nn_cnp, rn_fsp                   ! penalisation parameters 
    7784      !!---------------------------------------------------------------------- 
    7885      ! 
     
    97104      kpi = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells                 
    98105      kpj = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells               
    99  
    100 !!an 
    101 !#if defined key_agrif 
    102 !      IF( .NOT. Agrif_Root() ) THEN 
    103 !         kpi  = nbcellsx + 2 + 2*nbghostcells 
    104 !         kpj  = nbcellsy + 2 + 2*nbghostcells 
    105 !      ENDIF 
    106 !#endif 
    107 !!an 
    108106      ! 
    109107      IF( rn_modified_grav /= 0._wp) grav = rn_modified_grav   ! update the gravity  
    110       ! 
    111       !        !== Temporary namelist parameter ==! 
    112       ! 
    113       ! ll_dynldf_lap_PM = ln_dynldf_lap_PM 
    114108      ! 
    115109      kpk = jpkglo 
     
    117111      kperio = 0                    ! AM98 configuration : closed domain 
    118112      ! 
     113# if defined key_bvp 
     114      r1_abp = 1._wp / rn_abp 
     115#endif 
    119116      !                             ! control print 
    120117      IF(lwp) THEN 
     
    131128         WRITE(numout,*) '                    rotation angle chosen             rn_theta   = ', rn_theta, 'deg' 
    132129         WRITE(numout,*) '                    modified gravity          rn_modified_grav   = ', rn_modified_grav, 'm2/s' 
    133 #if defined key_agrif 
    134          IF( Agrif_Root() ) THEN 
    135 #endif 
    136          WRITE(numout,*) '         jpiglo = Nx*nn_AM98 + 2*n_ghost                    jpiglo = ', kpi 
    137          WRITE(numout,*) '         jpjglo = Nx*nn_AM98 + 2*n_ghost                    jpjglo = ', kpj 
    138 #if defined key_agrif 
    139          ENDIF 
    140 #endif 
    141130         WRITE(numout,*) '      number of model levels                              jpkglo = ', kpk 
    142131         WRITE(numout,*) '   ' 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90

    r13005 r13416  
    8888           ! length of the domain : 2000km x 2000km  
    8989           utau(ji,jj) = - ztauu * COS( rpi * gphiu(ji,jj) / 2000000._wp) 
    90            vtau(ji,jj) = - ztauv * COS( rpi * gphiu(ji,jj) / 2000000._wp) 
     90           vtau(ji,jj) = - ztauv * COS( rpi * gphiv(ji,jj) / 2000000._wp) 
    9191         END DO 
    9292      END DO 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90

    r13005 r13416  
    141141      pdepw_1d(1) =   0._wp 
    142142      pdept_1d(1) = 250._wp 
    143        
     143      ! 
    144144      pdepw_1d(2) = 500._wp 
    145145      pdept_1d(2) = 750._wp 
    146  
    147146      ! 
    148147      !                       ! e3t and e3w from depth 
     
    194193      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    195194      ! 
    196       ! rotated domain (45deg) 
    197       !zylim0 =  100000._wp * SQRT( 2._wp ) / 3    !  dy*sqrt(2)/2 * 2/3  
    198       !zylim1 = 2000000._wp - zylim0               !  2000km - dy*sqrt(2)/2 * 2/3 
    199       !zxlim0 =  100000._wp * SQRT( 2._wp ) / 3    !  dx*sqrt(2)/2 * 2/3  
    200       !zxlim1 = 2000000._wp - zxlim0               !  2000km - dx*sqrt(2)/2 * 2/3 
    201        
    202       ! rotated domain  
    203       !zylim0 =   10000._wp / REAL( nn_AM98, wp)      ! dy/10  
    204       !zylim1 = 2000000._wp + zylim0                  ! 2000km + dy/10  
    205       !zxlim0 =   10000._wp / REAL( nn_AM98, wp)      ! dx/10  
    206       !zxlim1 = 2000000._wp + zxlim0                  ! 2000km + dx/10 
    207  
    208195      !  
    209196      zylim0 =   10000._wp    ! +10km  
     
    211198      zxlim0 =   10000._wp    ! +10km  
    212199      zxlim1 = 2010000._wp    ! 2010km 
    213  
     200      ! 
    214201      DO jj = 1, jpj 
    215202         DO ji = 1, jpi 
     
    229216      END DO 
    230217      ! mask the lonely corners 
    231       DO jj = 1, jpj 
    232          DO ji = 1, jpi 
     218      DO jj = 2, jpim1 
     219         DO ji = 2, jpim1 
    233220         zcoeff = k_top(ji+1,jj) + k_top(ji,jj+1)   & 
    234221            +     k_top(ji-1,jj) + k_top(ji,jj-1) 
     
    239226         END DO 
    240227      END DO 
    241       ! k_bot(:,:) = NINT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
    242       ! 
    243       ! k_top(:,:) = MIN( 1 , k_bot(:,:) )      ! = 1    over the ocean point, =0 elsewhere 
     228      ! 
    244229      ! 
    245230   END SUBROUTINE zgr_msk_top_bot 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/cpp_AM98.fcm

    r12614 r13416  
    1  bld::tool::fppkeys  key_mpp_mpi key_iomput 
     1 bld::tool::fppkeys  key_mpp_mpi key_iomput key_qco 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/stpMLF.F90

    r13328 r13416  
    208208#endif 
    209209!!                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS 
    210                          CALL dyn_vor( kstp,     Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
     210                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    211211                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    212212      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90

    r13151 r13416  
    7575   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file 
    7676   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
     77   INTEGER ::   nid_F, nz_F, nh_F, ndim_F, ndim_hF   ! grid_F file 
    7778   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
    7879   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7980   INTEGER ::   ndex(1)                              ! ??? 
    80    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV, ndex_hF 
    8182   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    82    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
     83   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V, ndex_F 
    8384   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    8485 
     
    138139      CALL iom_put("e3u_0", e3u_0(:,:,:) ) 
    139140      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
     141      CALL iom_put("e3f_0", e3f_0(:,:,:) ) 
    140142      ! 
    141143!!st13 
     
    173175         CALL iom_put( "ssh" , ssh(:,:,Kmm) )                          ! sea surface height 
    174176      ENDIF 
    175  
    176 !!an  
    177         IF( iom_use("ht") )   &                  ! water column at t-point 
     177      ! 
     178      ! 
     179      IF( iom_use("ht") )   &                  ! water column at t-point 
    178180         CALL iom_put( "ht" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    179181      ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90

    r13151 r13416  
    1  
    21MODULE domvvl 
    32   !!====================================================================== 
     
    10421041               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 
    10431042         END_2D 
     1043!!an   dans le cas tourné, hf augmente et trend VOR diminue 
     1044!         DO_2D_10_10 
     1045!            zc3(ji,jj) =           (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  & 
     1046!               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  & 
     1047!               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  & 
     1048!               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) & 
     1049!               &           / MAX( tmask(ji,jj) + tmask(ji+1,jj) + tmask(ji,jj+1) + tmask(ji+1,jj+1), 1._wp )  
     1050!         END_2D 
     1051          
    10441052         CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 
    10451053         ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90

    r13151 r13416  
    2525   USE lib_mpp         ! distribued memory computing library 
    2626   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    27  
     27   ! 
     28   USE usrdef_nam , ONLY : ln_dynldf_lap_PM 
     29   ! 
    2830   IMPLICIT NONE 
    2931   PRIVATE 
     
    323325         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation  
    324326            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only) 
    325                ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
    326                WRITE(numout,*) ' ahmt tmask ' 
     327!!an          ! 
     328            WRITE(numout,*) '   ln_dynldf_lap_PM = ',ln_dynldf_lap_PM  
     329               IF(     ln_dynldf_lap_PM ) THEN                 !   laplacian operator (mask only) 
    327330!! mask tension at the coast (equivalent of ghostpoints at T) 
    328 !              DO jk = 1, jpk 
    329 !                 DO jj = 1, jpjm1 
    330 !                    DO ji = 1, jpim1      ! NO vector opt. 
    331 !                       ! si sum(fmask)==3 = mouillé (on touche pas) 
    332 !                       ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 
    333 !                       zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   & 
    334 !                          &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 
    335 !                       IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0 
    336 !                       ! 
    337 !                       !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   & 
    338 !                       !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 
    339 !                    END DO 
    340 !                 END DO 
    341 !              END DO 
    342 !              ahmt(jpi,:,1:jpkm1) =  0._wp 
    343 !              ahmt(:,jpj,1:jpkm1) =  0._wp 
    344 !              WRITE(numout,*) '   an45 ahmt x0' 
    345  
    346                ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
    347                WRITE(numout,*) ' ahmf fmask ' 
    348 !!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 
    349 !               DO jk = 1, jpkm1 
    350 !                  ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 
    351 !               END DO 
    352 !               WRITE(numout,*) '   an45 ahmf x2' 
    353  
     331                  DO jk = 1, jpk 
     332                     DO jj = 1, jpjm1 
     333                        DO ji = 1, jpim1      ! NO vector opt. 
     334                           ! si sum(fmask)==3 = mouillé (on touche pas) 
     335                           ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 
     336                           zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   & 
     337                              &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 
     338                           IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0 
     339                           ! 
     340                           !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   & 
     341                           !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 
     342                        END DO 
     343                     END DO 
     344                  END DO 
     345                  ahmt(jpi,:,1:jpkm1) =  0._wp 
     346                  ahmt(:,jpj,1:jpkm1) =  0._wp 
     347                  WRITE(numout,*) '  ahmt x0' 
     348!! apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 
     349                   DO jk = 1, jpkm1 
     350                      ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 
     351                   END DO 
     352                   WRITE(numout,*) '  ahmf x2' 
     353               ELSE 
     354               ! classic boundary condition on the viscosity coefficient 
     355                  ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
     356                  WRITE(numout,*) ' ahmt tmasked ' 
     357                  ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
     358                  WRITE(numout,*) ' ahmf fmasked ' 
     359               ENDIF 
     360!!an         !                  
    354361            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask) 
    355362               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 
Note: See TracChangeset for help on using the changeset viewer.