Changeset 13154 for NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90

Timestamp:

2020-06-24T15:31:32+02:00 (4 years ago)

Author:

gm

Message:

Alternative Directions on VOR and KE and bug fixed ln_vorlat

File:

• NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90 (modified) (3 diffs)

NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90

@@ -17 +17 @@
    USE phycst           ! physical constants
    USE usrdef_nam
-   !
+   USE lib_mpp        ! MPP library
    USE iom              ! xIOs server 
 
@@ -66 +66 @@
       REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
-      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
+      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva     
+      REAL(wp), DIMENSION(jpi,jpj)    ::   zssh
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3u, ze3v
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -141 +143 @@
                &            CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
-                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS
- 
-                            CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS
- 
-                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
 
 !!an - calcul du gradient de pression horizontal (explicit)
+      zssh(:,:) = ssh(:,:,Nnn)
+# if defined key_bvp
+      !   gradient and divergence not penalised
+      zssh(:,:) = zssh(:,:) * r1_rpo(:,:,1)  
+#endif 
       DO_3D_00_00( 1, jpkm1 )
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( zssh(ji+1,jj) - zssh(ji,jj) ) * r1_e1u(ji,jj)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( zssh(ji,jj+1) - zssh(ji,jj) ) * r1_e2v(ji,jj)
       END_3D
       !
-      ! add wind stress forcing and layer linear friction to the RHS 
+      
+!      IF( kstp == nit000 .AND. lwp ) THEN
+!         WRITE(numout,*)
+!         WRITE(numout,*) 'step.F90 : classic script used'
+!         WRITE(numout,*) '~~~~~~~'
+!         IF(       ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO   &
+!         &    .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO   ) THEN
+!            CALL ctl_stop('STOP','step : alternative direction asked but classis step'  )
+!         ENDIF
+!      ENDIF
+!!an     
+!                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)  ==> RHS
+! 
+!                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity             ==> RHS
+! 
+!!an     In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90
+  
+                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs)  ! vorticity            ==> RHS
+  
+                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
+
+      ! add wind stress forcing and layer linear friction to the RHS
+      ze3u(:,:,:) = e3u(:,:,:,Nnn)
+      ze3v(:,:,:) = e3v(:,:,:,Nnn)
+# if defined key_bvp
+      !ze3u(:,:,:) = ze3u(:,:,:) * r1_rpo(:,:,:)  
+      !ze3v(:,:,:) = ze3v(:,:,:) * r1_rpo(:,:,:)  
+#endif  
       z1_2rho0 = 0.5_wp * r1_rho0
       DO_3D_00_00(1,jpkm1)
-         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / ze3u(ji,jj,jk)   &
             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
-            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ze3v(ji,jj,jk)   &
+            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)  
       END_3D   
-!!an         
+      !
+# if defined key_bvp
+      !  Add frictionnal term   - sigma * u
+      ! can be done via sigmaT (mean_ij)
+      DO_3D_00_00(1,jpkm1)
+         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - bmpu(ji,jj,jk) * uu(ji,jj,jk,Nbb)
+         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - bmpv(ji,jj,jk) * vv(ji,jj,jk,Nbb)
+      END_3D 
+#endif 
+      !         
   
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>