Changeset 11394 for branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90

Timestamp:

2019-08-02T15:14:02+02:00 (5 years ago)

Author:

mattmartin

Message:

First implementation of STOPACK in the GO6 package branch.

File:

• branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90 (modified) (13 diffs)

branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -118 +118 @@
          CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. )
       END DO
-
+      
+      IF( ln_stopack )   CALL stopack_pert( kstp )
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
+
+      IF( ln_sto_eos )   CALL sto_par( kstp )          ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( tsn  )          ! Random T/S fluctuations
+      IF( ln_skeb  )     CALL skeb_comp( kstp )
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -149 +153 @@
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
+         IF ( nn_spp_arnf .GT. 0 ) THEN
+              rn_avt_rnf0 = rn_avt_rnf
+              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf )
+         ENDIF
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
       IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
@@ -163 +171 @@
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
       IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
+      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' )
       !
       !  LATERAL  PHYSICS
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+         IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
          IF( ln_zps .AND. .NOT. ln_isfcav)                               &
             &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
@@ -202 +212 @@
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
-            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
+                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                               &
                &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
@@ -214 +225 @@
                                   va(:,:,:) = 0.e0
           IF(  lk_asminc .AND. ln_asmiau .AND. &
-             & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified)
-          IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends
-                                  CALL dyn_adv      ( kstp )   ! advection (vector or flux form)
-                                  CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis
-                                  CALL dyn_ldf      ( kstp )   ! lateral mixing
-          IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified)
+             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment
+          IF( ln_sppt_dyn )       CALL dyn_sppt_apply( kstp, 0 )
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified)
+          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends
+                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form)
+                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis
+                                  CALL dyn_ldf       ( kstp )   ! lateral mixing
+          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified)
 #if defined key_agrif
           IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge
@@ -262 +274 @@
                              tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero
 
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
+      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 0 )
+      IF( lk_asminc .AND. ln_asmiau .AND. &
+        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -280 +293 @@
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
 #endif
+      IF( ln_sppt_tra )      CALL tra_sppt_apply ( kstp, 1 )
                              CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields
 
@@ -285 +299 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+         IF( ln_sppt_tra )   CALL tra_sppt_apply ( kstp, 2 )
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
             IF( ln_zps .AND. .NOT. ln_isfcav)                                &
@@ -308 +323 @@
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
          IF( ln_bias )       CALL dyn_bias( kstp )
+         IF( ln_sppt_tra )   THEN
+                                CALL tra_sppt_apply ( kstp, 2 )
+                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
+                                IF( ln_zps .AND. .NOT. ln_isfcav)                &
+                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+                                IF( ln_zps .AND.       ln_isfcav)                                   & 
+                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+                                &                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
+         ENDIF
       ENDIF
 
@@ -324 +350 @@
 
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
       ELSE
@@ -331 +358 @@
         IF(  lk_asminc .AND. ln_asmiau .AND. &
            & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+        IF( ln_sppt_dyn )      CALL dyn_sppt_apply ( kstp, 0 )
         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
@@ -343 +371 @@
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                CALL dyn_bfr( kstp )         ! bottom friction
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 1 )
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
       ENDIF
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
+            IF( ln_sppt_dyn )  CALL dyn_sppt_apply ( kstp, 2 )
+            IF( ln_skeb )      CALL skeb_apply ( kstp )
 
                                CALL ssh_swp( kstp )         ! swap of sea surface height
@@ -381 +412 @@
          IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
       ENDIF
+      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters