Changeset 13088

Timestamp:

2020-06-10T13:13:39+02:00 (5 years ago)

Author:

jwhile

Message:

Bug fixes for 1D running

Location:

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC

Files:

• ASM/asmbgc.F90 (modified) (53 diffs)
• DYN/sshwzv.F90 (modified) (15 diffs)
• SBC/albedo.F90 (modified) (18 diffs)
• SBC/sbcblk_core.F90 (modified) (34 diffs)
• SBC/sbcssr.F90 (modified) (8 diffs)
• STO/stopack.F90 (modified) (47 diffs)
• TRA/trabbc.F90 (modified) (6 diffs)
• TRA/trabbl.F90 (modified) (6 diffs)
• TRA/traqsr.F90 (modified) (33 diffs)
• ZDF/zdfbfr.F90 (modified) (4 diffs)
• ZDF/zdfevd.F90 (modified) (7 diffs)
• ZDF/zdfgls.F90 (modified) (39 diffs)
• ZDF/zdftke.F90 (modified) (37 diffs)
• step.F90 (modified) (10 diffs)

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

@@ -33 +33 @@
    USE oce, ONLY:           & ! active tracer variables
       & tsn
-   USE zdfmxl, ONLY :       & ! mixed layer depth 
+   USE zdfmxl, ONLY :       & ! mixed layer depth
 #if defined key_karaml
       & hmld_kara,          &
       & ln_kara,            &
-#endif   
-      & hmld,               & 
+#endif
+      & hmld,               &
       & hmlp,               &
-      & hmlpt 
+      & hmlpt
    USE asmpar, ONLY:        & ! assimilation parameters
       & c_asmbkg,           &
@@ -89 +89 @@
 
    IMPLICIT NONE
-   PRIVATE                   
+   PRIVATE
 
    PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90
@@ -290 +290 @@
 
       ! Allocate and read increments
-      
+
       IF ( ln_slchltotinc ) THEN
          ALLOCATE( slchltot_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
       ENDIF
-      
+
       IF ( ln_slchldiainc ) THEN
          ALLOCATE( slchldia_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
       ENDIF
-      
+
       IF ( ln_slchlnoninc ) THEN
          ALLOCATE( slchlnon_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
       ENDIF
-      
+
       IF ( ln_schltotinc ) THEN
          ALLOCATE( schltot_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
       ENDIF
-      
+
       IF ( ln_slphytotinc ) THEN
          ALLOCATE( slphytot_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
       ENDIF
-      
+
       IF ( ln_slphydiainc ) THEN
          ALLOCATE( slphydia_bkginc(jpi,jpj) )
          CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
       ENDIF
-      
+
       IF ( ln_slphynoninc ) THEN
          ALLOCATE( slphynon_bkginc(jpi,jpj) )
@@ -335 +335 @@
          CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
       ENDIF
-      
+
       IF ( ln_plchltotinc ) THEN
          ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
       ENDIF
-      
+
       IF ( ln_pchltotinc ) THEN
          ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
       ENDIF
-      
+
       IF ( ln_pno3inc ) THEN
          ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
       ENDIF
-      
+
       IF ( ln_psi4inc ) THEN
          ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
       ENDIF
-      
+
       IF ( ln_pdicinc ) THEN
          ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
       ENDIF
-      
+
       IF ( ln_palkinc ) THEN
          ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
       ENDIF
-      
+
       IF ( ln_pphinc ) THEN
          ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
          CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
       ENDIF
-      
+
       IF ( ln_po2inc ) THEN
          ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
@@ -377 +377 @@
 
       ! Allocate balancing increments
-      
+
       IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
@@ -388 +388 @@
 #endif
       ENDIF
-      
+
       IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
 #if defined key_top
@@ -443 +443 @@
       ! Initialise
       p_incs(:,:) = 0.0
-      
+
       ! read from file
       CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
-      
+
       ! Apply the masks
       p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
-      
+
       ! Set missing increments to 0.0 rather than 1e+20
       ! to allow for differences in masks
@@ -481 +481 @@
       ! Initialise
       p_incs(:,:,:) = 0.0
-      
+
       ! read from file
       CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
-      
+
       ! Apply the masks
       p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
-      
+
       ! Set missing increments to 0.0 rather than 1e+20
       ! to allow for differences in masks
@@ -538 +538 @@
          cchl_p_bkg(:,:,:) = 0.0
 #endif
-         
+
          !--------------------------------------------------------------------
          ! Read background variables for phytoplankton assimilation
@@ -558 +558 @@
          CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
 #endif
-         
+
          IF ( ln_phytobal ) THEN
 
@@ -602 +602 @@
 
          CALL iom_close( inum )
-         
+
          DO jt = 1, jptra
             tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
          END DO
-      
+
       ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
 
@@ -615 +615 @@
 
          CALL iom_open( c_asmbkg, inum )
-         
+
 #if defined key_hadocc
          CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
@@ -626 +626 @@
 
          CALL iom_close( inum )
-         
+
          DO jt = 1, jptra
             tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
          END DO
          mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
-      
+
       ENDIF
 #else
@@ -655 +655 @@
       !!
       !! ** Action  :
-      !!                   
+      !!
       !! References :
       !!
@@ -669 +669 @@
       REAL(wp)          :: zdate         ! Date
       !!------------------------------------------------------------------------
-     
+
       ! Set things up
       zdate = REAL( ndastp )
@@ -680 +680 @@
             &                    TRIM( c_asmbal ) // ' at timestep = ', kt
 
-         ! Define the output file       
+         ! Define the output file
          CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
 
@@ -767 +767 @@
             &                    TRIM( c_asmbal ) // ' at timestep = ', kt
       ENDIF
-                                 
+
    END SUBROUTINE asm_bgc_bal_wri
 
@@ -847 +847 @@
       !! ** Action  :   return non-log increments
       !!
-      !! References :   
+      !! References :
       !!------------------------------------------------------------------------
       !!
@@ -879 +879 @@
       !!------------------------------------------------------------------------
       !!                    ***  ROUTINE phyto2d_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the chlorophyll assimilation increments.
       !!
@@ -886 +886 @@
       !!              Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!------------------------------------------------------------------------
       INTEGER,  INTENT(IN) :: kt        ! Current time step
@@ -914 +914 @@
 #endif
       !!------------------------------------------------------------------------
-      
+
       IF ( kt <= nit000 ) THEN
 
@@ -920 +920 @@
          ! Remember that two sets of non-log increments should not be
          ! expected to be in the same ratio as their log equivalents
-         
+
          ! Total chlorophyll
          IF ( ln_slchltotinc ) THEN
@@ -1074 +1074 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
                   &  kt,' with IAU weight = ', pwgtiau(it)
@@ -1105 +1105 @@
          ENDIF
 
-      ELSEIF ( ll_asmdin ) THEN 
+      ELSEIF ( ll_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -1134 +1134 @@
             END WHERE
 #endif
- 
+
             ! Do not deallocate arrays - needed by asm_bgc_bal_wri
             ! which is called at end of model run
@@ -1150 +1150 @@
       !!------------------------------------------------------------------------
       !!                    ***  ROUTINE phyto3d_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the profile chlorophyll assimilation increments.
       !!
@@ -1156 +1156 @@
       !!              Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!------------------------------------------------------------------------
       INTEGER,  INTENT(IN) :: kt        ! Current time step
@@ -1261 +1261 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
                   &  kt,' with IAU weight = ', pwgtiau(it)
@@ -1292 +1292 @@
          ENDIF
 
-      ELSEIF ( ll_asmdin ) THEN 
+      ELSEIF ( ll_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -1321 +1321 @@
             END WHERE
 #endif
- 
+
             ! Do not deallocate arrays - needed by asm_bgc_bal_wri
             ! which is called at end of model run
@@ -1338 +1338 @@
       !!------------------------------------------------------------------------
       !!                    ***  ROUTINE pco2_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the pco2/fco2 assimilation increments.
       !!
@@ -1345 +1345 @@
       !!              Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!------------------------------------------------------------------------
       INTEGER, INTENT(IN)                   :: kt           ! Current time step
@@ -1495 +1495 @@
                jkmax = jpk-1
                DO jk = jpk-1, 1, -1
+#if defined key_vvl
                   IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
                      & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
@@ -1500 +1501 @@
                      jkmax = jk
                   ENDIF
+#else
+                  IF ( ( zmld(ji,jj) >  gdepw_0(ji,jj,jk)   ) .AND. &
+                     & ( zmld(ji,jj) <= gdepw_0(ji,jj,jk+1) ) ) THEN
+                     zmld(ji,jj) = gdepw_0(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+#endif
                END DO
                !
@@ -1526 +1534 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                IF ( ln_spco2inc ) THEN
                   WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
@@ -1563 +1571 @@
          ENDIF
 
-      ELSEIF ( ll_asmdin ) THEN 
+      ELSEIF ( ll_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -1592 +1600 @@
             END WHERE
 #endif
- 
+
             ! Do not deallocate arrays - needed by asm_bgc_bal_wri
             ! which is called at end of model run
@@ -1609 +1617 @@
       !!------------------------------------------------------------------------
       !!                    ***  ROUTINE ph_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the pH assimilation increments.
       !!
@@ -1616 +1624 @@
       !!              Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!------------------------------------------------------------------------
       INTEGER,                          INTENT(IN) :: kt        ! Current time step
@@ -1626 +1634 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
-      
+
       REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
       REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
@@ -1679 +1687 @@
             sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
          ENDIF
-         
+
          ! Account for pCO2 balancing if required
          IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
@@ -1685 +1693 @@
             alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
          ENDIF
-         
+
          ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
          ! This requires three calls to the MOCSY carbonate package
@@ -1750 +1758 @@
                      ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
                      ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
-                     
+
                   ENDIF
-                  
+
                END DO
             END DO
@@ -1758 +1766 @@
 
       ENDIF
-      
+
       IF ( ll_asmiau ) THEN
 
@@ -1772 +1780 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
                   &  kt,' with IAU weight = ', pwgtiau(it)
@@ -1794 +1802 @@
          ENDIF
 
-      ELSEIF ( ll_asmdin ) THEN 
+      ELSEIF ( ll_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -1814 +1822 @@
                trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
             END WHERE
- 
+
             ! Do not deallocate arrays - needed by asm_bgc_bal_wri
             ! which is called at end of model run
@@ -1820 +1828 @@
          !
       ENDIF
-#endif      
+#endif
       !
    END SUBROUTINE ph_asm_inc
@@ -1831 +1839 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dyn_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER,  INTENT(IN) :: kt        ! Current time step
@@ -1983 +1991 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
                   &  kt,' with IAU weight = ', pwgtiau(it)
@@ -2071 +2079 @@
 #endif
             ENDIF
-           
+
             IF ( kt == nitiaufin_r ) THEN
                IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
@@ -2082 +2090 @@
          ENDIF
 
-      ELSEIF ( ll_asmdin ) THEN 
+      ELSEIF ( ll_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -2179 +2187 @@
 #endif
             ENDIF
- 
+
             IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
             IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -1 @@
-MODULE sshwzv   
+MODULE sshwzv
    !!==============================================================================
    !!                       ***  MODULE  sshwzv  ***
@@ -5 +5 @@
    !!==============================================================================
    !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
    !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
@@ -17 +17 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
+   USE dom_oce         ! ocean space and time domain variables
    USE sbc_oce         ! surface boundary condition: ocean
    USE domvvl          ! Variable volume
@@ -28 +28 @@
    USE lib_mpp         ! MPP library
    USE bdy_oce
-   USE bdy_par         
+   USE bdy_par
    USE bdydyn2d        ! bdy_ssh routine
 #if defined key_agrif
    USE agrif_opa_interp
 #endif
-#if defined key_asminc   
+#if defined key_asminc
    USE asminc          ! Assimilation increment
 #endif
@@ -59 +59 @@
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE ssh_nxt  ***
-      !!                   
+      !!
       !! ** Purpose :   compute the after ssh (ssha)
       !!
@@ -73 +73 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
       INTEGER, INTENT(in) ::   kt                      ! time step
-      ! 
+      !
       INTEGER             ::   jk                      ! dummy loop indices
       REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
@@ -80 +80 @@
       IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
       !
-      CALL wrk_alloc( jpi, jpj, zhdiv ) 
+      CALL wrk_alloc( jpi, jpj, zhdiv )
       !
       IF( kt == nit000 ) THEN
@@ -102 +102 @@
 #if defined key_vvl
 ! Don't directly adjust ssh but change hdivn at all levels instead
-! In trasbc also add in the heat and salt content associated with these changes at each level  
-        DO jk = 1, jpkm1                                 
-                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
+! In trasbc also add in the heat and salt content associated with these changes at each level
+        DO jk = 1, jpkm1
+                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk)
         END DO
-      ENDIF
-#endif
+#endif
+      ENDIF
 #endif
 
@@ -121 +121 @@
       ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
-      ! 
+      !
       z1_rau0 = 0.5_wp * r1_rau0
       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
@@ -146 +146 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
       !
-      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
+      CALL wrk_dealloc( jpi, jpj, zhdiv )
       !
       IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt')
@@ -152 +152 @@
    END SUBROUTINE ssh_nxt
 
-   
+
    SUBROUTINE wzv( kt )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wzv  ***
-      !!                   
+      !!
       !! ** Purpose :   compute the now vertical velocity
       !!
-      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
-      !!      velocity is computed by integrating the horizontal divergence  
+      !! ** Method  : - Using the incompressibility hypothesis, the vertical
+      !!      velocity is computed by integrating the horizontal divergence
       !!      from the bottom to the surface minus the scale factor evolution.
       !!        The boundary conditions are w=0 at the bottom (no flux) and.
@@ -176 +176 @@
       REAL(wp)            ::   z1_2dt       ! local scalars
       !!----------------------------------------------------------------------
-      
+
       IF( nn_timing == 1 )  CALL timing_start('wzv')
       !
@@ -195 +195 @@
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
-         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
+         CALL wrk_alloc( jpi, jpj, jpk, zhdiv )
          !
          DO jk = 1, jpkm1
@@ -215 +215 @@
          END DO
          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
-         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
+         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )
       ELSE   ! z_star and linear free surface cases
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
@@ -241 +241 @@
       !!                    ***  ROUTINE ssh_nxt  ***
       !!
-      !! ** Purpose :   achieve the sea surface  height time stepping by 
+      !! ** Purpose :   achieve the sea surface  height time stepping by
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_nxt  
+      !!              ssha  already computed in ssh_nxt
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

@@ -21 +21 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE stopack
 
@@ -31 +31 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-  
+
    REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
    REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
@@ -37 +37 @@
    REAL(wp) ::   c2       = 0.10    !  "        "
    REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
- 
+
    !                             !!* namelist namsbc_alb
    INTEGER  ::   nn_ice_alb
@@ -52 +52 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_ice  ***
-      !!          
-      !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!       
+      !!
+      !! ** Purpose :   Computation of the albedo of the snow/ice system
+      !!
       !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
       !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
@@ -73 +73 @@
       !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
       !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
-      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger
       !!                     under melting conditions than under freezing conditions
-      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows
       !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
       !!
       !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
       !!                Brandt et al. 2005, J. Climate, vol 18
-      !!                Grenfell & Perovich 2004, JGR, vol 109 
+      !!                Grenfell & Perovich 2004, JGR, vol 109
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -97 +97 @@
 
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-      
+
       CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
-      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
-
-      
+      IF( albd_init == 0 )   CALL albedo_init      ! initialization
+
+
       SELECT CASE ( nn_ice_alb )
 
@@ -109 +109 @@
       !------------------------------------------
       CASE( 0 )
-       
+
          ralb_sf = 0.80       ! dry snow
          ralb_sm = 0.65       ! melting snow
          ralb_if = 0.72       ! bare frozen ice
-         ralb_im = rn_albice  ! bare puddled ice 
-         
+         ralb_im = rn_albice  ! bare puddled ice
+
          !  Computation of ice albedo (free of snow)
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
          ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
          END  WHERE
-      
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
@@ -127 +127 @@
          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
-     
+
          DO jl = 1, ijpl
             DO jj = 1, jpj
@@ -133 +133 @@
                   ! freezing snow
                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  !                                        !  freezing snow        
+                  !                                        !  freezing snow
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
                   zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
                      &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
-                     &        +         zswitch   * ralb_sf  
+                     &        +         zswitch   * ralb_sf
 
                   ! melting snow
@@ -143 +143 @@
                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
                   zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
-                      &     +         zswitch   *   ralb_sm 
+                      &     +         zswitch   *   ralb_sm
                   !
                   ! snow albedo
-                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
-               
+
                   ! Ice/snow albedo
                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
@@ -155 +155 @@
                END DO
             END DO
-            
+
+#if defined key_traldf_c2d || key_traldf_c3d
             IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
                & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
-                        
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
 
@@ -166 +171 @@
       !  New parameterization (2016)
       !------------------------------------------
-      CASE( 1 ) 
+      CASE( 1 )
 
          ralb_im = rn_albice  ! bare puddled ice
@@ -181 +186 @@
 !         ralb_sm = 0.82      ! melting snow
 !         ralb_if = 0.54      ! bare frozen ice
-! 
+!
          !  Computation of ice albedo (free of snow)
-         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )
          z1_c2 = 1. / 0.05
          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
          ELSE WHERE                                              ;   zalb = ralb_if
          END  WHERE
-         
+
          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
@@ -205 +210 @@
 
                    ! snow albedo
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )
                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
 
-                  ! Ice/snow albedo   
+                  ! Ice/snow albedo
                   zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
                   pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
@@ -214 +219 @@
               END DO
             END DO
-            
+
+#if defined key_traldf_c2d || key_traldf_c3d
             IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
                & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl )
-            
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
          END DO
          ! Effect of the clouds (2d order polynomial)
-         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );
 
       END SELECT
-      
+
       CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
@@ -232 +242 @@
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE albedo_oce  ***
-      !! 
+      !!
       !! ** Purpose :   Computation of the albedo of the ocean
       !!----------------------------------------------------------------------
@@ -238 +248 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) :: zcoef 
+      REAL(wp) :: zcoef
       !!----------------------------------------------------------------------
       !
       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
-      pa_oce_cs(:,:) = zcoef 
+      pa_oce_cs(:,:) = zcoef
       pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
@@ -257 +267 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice
       !!----------------------------------------------------------------------
       !

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -161 +161 @@
          !                                      ! ====================== !
          !
-         rn_sfac = 1._wp       ! Default to one if missing from namelist 
+         rn_sfac = 1._wp       ! Default to one if missing from namelist
          REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
          READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
@@ -205 +205 @@
       ENDIF
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN
          rn_vfac0(:,:) = rn_vfac
          CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw )
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_relw > 0 ) &
+         & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
@@ -217 +223 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
+         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)
          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
          ENDIF
-         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
+         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
@@ -231 +237 @@
       !
    END SUBROUTINE sbc_blk_core
-   
-   
+
+
    SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv )
       !!---------------------------------------------------------------------
@@ -242 +248 @@
       !! ** Method  :   CORE bulk formulea for the ocean using atmospheric
       !!      fields read in sbc_read
-      !! 
+      !!
       !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
       !!              - vtau    : j-component of the stress at V-point  (N/m2)
@@ -280 +286 @@
       ! local scalars ( place there for vector optimisation purposes)
       zcoef_qsatw = 0.98 * 640380. / rhoa
-      
+
       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
 
@@ -288 +294 @@
 
       ! ... components ( U10m - U_oce ) at T-point (unmasked)
-      zwnd_i(:,:) = 0.e0  
+      zwnd_i(:,:) = 0.e0
       zwnd_j(:,:) = 0.e0
 #if defined key_cyclone
@@ -332 +338 @@
       CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   &
          &               Cd, Ch, Ce, zt_zu, zq_zu )
-    
+
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -363 +369 @@
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
 
-    
+
       !  Turbulent fluxes over ocean
       ! -----------------------------
@@ -388 +394 @@
          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ')
       ENDIF
-       
+
       ! ----------------------------------------------------------------------------- !
       !     III    Total FLUXES                                                       !
@@ -396 +402 @@
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
       !
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
@@ -437 +443 @@
       !
    END SUBROUTINE blk_oce_core
- 
-   
+
+
 #if defined key_lim2 || defined key_lim3
    SUBROUTINE blk_ice_core_tau
@@ -531 +537 @@
 
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
-      
+
    END SUBROUTINE blk_ice_core_tau
 
@@ -544 +550 @@
       !!                between atmosphere and sea-ice using CORE bulk
       !!                formulea, ice variables and read atmmospheric fields.
-      !! 
+      !!
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
@@ -564 +570 @@
       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx')
       !
-      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
+      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
 
       ! local scalars ( place there for vector optimisation purposes)
@@ -587 +593 @@
                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
-               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
+               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
 
                ! ----------------------------!
@@ -597 +603 @@
                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
                ! Latent Heat
-               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &
                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
               ! Latent heat sensitivity for ice (Dqla/Dt)
@@ -628 +634 @@
 
 #if defined  key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_alloc( jpi,jpj, zevap, zsnw )
 
       ! --- evaporation --- !
@@ -638 +644 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
+      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
@@ -661 +667 @@
       DO jl = 1, jpl
          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
-                                   ! But we do not have Tice => consider it at 0 degC => evap=0 
+                                   ! But we do not have Tice => consider it at 0 degC => evap=0
       END DO
 
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
+      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )
 #endif
 
@@ -688 +694 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
-      
+
    END SUBROUTINE blk_ice_core_flx
 #endif
@@ -701 +707 @@
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
       !!
-      !! ** Method : Monin Obukhov Similarity Theory 
+      !! ** Method : Monin Obukhov Similarity Theory
       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
       !!
@@ -711 +717 @@
       !!    - better first guess of stability by checking air-sea difference of virtual temperature
       !!       rather than temperature difference only...
-      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
+      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
@@ -746 +752 @@
 
       IF( nn_timing == 1 )  CALL timing_start('turb_core_2z')
-    
+
       CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd )
       CALL wrk_alloc( jpi,jpj, zeta_u, stab )
@@ -758 +764 @@
       U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
 
-      !! First guess of stability: 
+      !! First guess of stability:
       ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt
       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
@@ -772 +778 @@
       Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 )
       Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
-    
+
       !! Initializing transf. coeff. with their first guess neutral equivalents :
       Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10
@@ -783 +789 @@
          !
          ztmp1 = T_zu - sst   ! Updating air/sea differences
-         ztmp2 = q_zu - q_sat 
+         ztmp2 = q_zu - q_sat
 
          ! Updating turbulent scales :   (L&Y 2004 eq. (7))
          ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta*
          ztmp2  = Ce/sqrt_Cd*ztmp2    ! q*
-       
+
          ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu
 
          ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
-         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
+         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)
          !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu)
 
@@ -799 +805 @@
          zpsi_h_u = psi_h( zeta_u )
          zpsi_m_u = psi_m( zeta_u )
-       
+
          !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
          IF ( .NOT. l_zt_equal_zu ) THEN
@@ -808 +814 @@
             q_zu = max(0., q_zu)
          END IF
-       
+
          IF( ln_cdgw ) THEN      ! surface wave case
-            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
+            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )
             Cd      = sqrt_Cd * sqrt_Cd
          ELSE
@@ -820 +826 @@
            ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10
            sqrt_Cd_n10 = sqrt(ztmp0)
-       
+
            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b)
            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
@@ -827 +833 @@
            !! Update of transfer coefficients:
            ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a)
-           Cd      = ztmp0 / ( ztmp1*ztmp1 )   
+           Cd      = ztmp0 / ( ztmp1*ztmp1 )
            sqrt_Cd = SQRT( Cd )
          ENDIF
@@ -833 +839 @@
          ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
          ztmp2 = sqrt_Cd / sqrt_Cd_n10
-         ztmp1 = 1. + Ch_n10*ztmp0               
+         ztmp1 = 1. + Ch_n10*ztmp0
          Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
          !
-         ztmp1 = 1. + Ce_n10*ztmp0               
+         ztmp1 = 1. + Ce_n10*ztmp0
          Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
          !
@@ -854 +860 @@
    FUNCTION cd_neutral_10m( zw10 )
       !!----------------------------------------------------------------------
-      !! Estimate of the neutral drag coefficient at 10m as a function 
+      !! Estimate of the neutral drag coefficient at 10m as a function
       !! of neutral wind  speed at 10m
       !!
@@ -860 +866 @@
       !!
       !! Author: L. Brodeau, june 2014
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s)
       REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m
       !
       REAL(wp), DIMENSION(:,:), POINTER ::   rgt33
-      !!----------------------------------------------------------------------    
+      !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj, rgt33 )
       !
       !! When wind speed > 33 m/s => Cyclone conditions => special treatment
-      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1  
+      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1
       cd_neutral_10m = 1.e-3 * ( &
          &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33.

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

@@ -24 +24 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE stopack
    USE wrk_nemo       ! Memory Allocation
@@ -42 +42 @@
    REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS
    REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS
-   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
+   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
 
@@ -101 +101 @@
                rn_dqdt_s(:,:) = rn_dqdt
 
-               IF( ln_stopack .AND. nn_spp_dqdt > 0 ) &
-                  & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt )
+#if defined key_traldf_c2d || key_traldf_c3d
+            IF( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt )
+#else
+            IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
                DO jj = 1, jpj
                   DO ji = 1, jpi
@@ -117 +124 @@
                CALL wrk_alloc( jpi, jpj, zsrp)
                zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
                IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
                   & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
+#else
+            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) &
+               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
 !CDIR COLLAPSE
                DO jj = 1, jpj
                   DO ji = 1, jpi
                      zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
-                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
+                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
                      erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only)
@@ -134 +149 @@
                CALL wrk_alloc( jpi, jpj, zsrp)
                zsrp(:,:) = rn_deds
+#if defined key_traldf_c2d || key_traldf_c3d
                IF( ln_stopack .AND. nn_spp_dedt > 0 ) &
                   & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds )
-               zerp_bnd = rn_sssr_bnd / rday                          !       -              -    
+#else
+               IF ( ln_stopack .AND. nn_spp_dedt > 0 ) &
+                  & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// &
+                                   'key_traldf_c2d or key_traldf_c3d')
+#endif
+               zerp_bnd = rn_sssr_bnd / rday                          !       -              -
 !CDIR COLLAPSE
                DO jj = 1, jpj
-                  DO ji = 1, jpi                            
+                  DO ji = 1, jpi
                      zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
@@ -161 +182 @@
    END SUBROUTINE sbc_ssr
 
- 
+
    SUBROUTINE sbc_ssr_init
       !!---------------------------------------------------------------------
@@ -184 +205 @@
       !!----------------------------------------------------------------------
       !
- 
-      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist : 
+
+      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :
       READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp )
@@ -240 +261 @@
       ENDIF
       !
-      !                            !* Initialize qrp and erp if no restoring 
+      !                            !* Initialize qrp and erp if no restoring
       IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
       IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
       !
    END SUBROUTINE sbc_ssr_init
-      
+
    !!======================================================================
 END MODULE sbcssr

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90

@@ -46 +46 @@
 #define numnam_ref numnam
 #define numnam_cfg numnam
-#define lwm        lwp   
+#define lwm        lwp
 #define numond     numout
 
-#define wmask      tmask  
+#define wmask      tmask
 
 #endif
@@ -63 +63 @@
    !!                                          (SPP, SKEB and sea-ice)
    !!----------------------------------------------------------------------
-   !!   
+   !!
    !!   stopack       : Generate stochastic physics perturbations
-   !!   
+   !!
    !!                   Method
    !!                   ======
@@ -71 +71 @@
    !!       - SPPT (Stochastically perturbed parameterization
    !!         tendencies )scheme for user-selected trends for
-   !!            tracers, momentum and sea-ice 
+   !!            tracers, momentum and sea-ice
    !!       - SPP (Schastically perturbed parameters) scheme
    !!         for some (namelist) parameters
@@ -77 +77 @@
    !!         backscatter energy dissipated numerically or
    !!            through deep convection.
-   !!   
-   !!   
+   !!
+   !!
    !!                   Acknowledgements: C3S funded ERGO project
-   !!   
+   !!
    !!----------------------------------------------------------------------
    USE par_kind
    USE timing          ! Timing
    USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
+   USE dom_oce         ! ocean space and time domain variables
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -104 +104 @@
    USE wrk_nemo
    USE diaptr
-   USE zdf_oce         
+   USE zdf_oce
    USE phycst
 
@@ -113 +113 @@
    PUBLIC tra_sppt_collect
    PUBLIC dyn_sppt_collect
-   PUBLIC tra_sppt_apply 
-   PUBLIC dyn_sppt_apply 
+   PUBLIC tra_sppt_apply
+   PUBLIC dyn_sppt_apply
    PUBLIC stopack_rst
    PUBLIC stopack_init
@@ -140 +140 @@
    REAL(wp), SAVE :: rn_spp_tau, rn_spp_stdev
    INTEGER :: skeb_filter_pass, spp_filter_pass
- 
+
    ! SPPT Logical switches for individual tendencies
    LOGICAL :: ln_sppt_taumap, ln_stopack_restart, ln_distcoast, &
    ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
    ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
-   ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf    
+   ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf
    LOGICAL :: &
    ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad,&
@@ -181 +181 @@
    INTEGER, PARAMETER, PUBLIC :: jk_spp_qsi0   = 8
    INTEGER, PARAMETER, PUBLIC :: jk_spp_bfr    = 9
-   INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd   = 10 
+   INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd   = 10
    INTEGER, PARAMETER, PUBLIC :: jk_spp_avt    = 11
    INTEGER, PARAMETER, PUBLIC :: jk_spp_avm    = 12
@@ -219 +219 @@
    INTEGER, SAVE :: numrep        = 602
    INTEGER, SAVE :: lkt
-   
+
    ! Randome generator seed
    INTEGER, SAVE   :: nn_stopack_seed(4)
@@ -275 +275 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $ 
+   !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -289 +289 @@
       !!
       !! ** Purpose :   Collect tracer tendencies (additive)
-      !!                This function is called by the tendency diagnostics 
+      !!                This function is called by the tendency diagnostics
       !!                module
       !!----------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
 #endif
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
@@ -354 +354 @@
       !!
       !! ** Purpose :   Collect momentum tendencies (additive)
-      !!                This function is called by the tendency diagnostics 
+      !!                This function is called by the tendency diagnostics
       !!                module
       !!----------------------------------------------------------------------
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
 #endif
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity
       INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
@@ -469 +469 @@
       !!
       !! ** Purpose :   Apply collinear perturbation to ice fields
-      !!                For specific processes coded in LIM2/LIM3 
+      !!                For specific processes coded in LIM2/LIM3
       !!----------------------------------------------------------------------
       !
@@ -529 +529 @@
          zicewrk(:,:,jm) = z5     ; jm=jm+1
          zicewrk(:,:,jm) = z6     ; jm=jm+1
-         zicewrk(:,:,jm) = z7   
+         zicewrk(:,:,jm) = z7
        ENDIF
        IF( kopt .EQ. 3 ) THEN
@@ -601 +601 @@
             CALL lbc_lnk( u_ice, 'U', -1. )
             CALL lbc_lnk( v_ice, 'V', -1. )
-#endif         
+#endif
 #if defined key_lim2   &&   defined key_lim2_vp
             CALL lbc_lnk( u_ice(:,1:jpj), 'I', -1. )
             CALL lbc_lnk( v_ice(:,1:jpj), 'I', -1. )
-#endif         
+#endif
           ENDIF
           DEALLOCATE ( zicewrk )
@@ -616 +616 @@
       !!                  ***  ROUTINE spp_gen ***
       !!
-      !! ** Purpose :   Perturbing parameters (generic function)  
+      !! ** Purpose :   Perturbing parameters (generic function)
       !!                Given a value of standard deviation, the 2D parameter
       !!                coeff is perturbed following an additive Normal,
@@ -624 +624 @@
       !!----------------------------------------------------------------------
    INTEGER, INTENT( in ) :: kt
+#if defined key_traldf_c3d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
+   REAL(wp), POINTER, DIMENSION(:,:,:) ::   gauss
+#elif defined key_traldf_c2d
    REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
+   REAL(wp), POINTER, DIMENSION(:,:) ::   gauss
+#elif defined key_traldf_c1d
+   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
+   REAL(wp), POINTER, DIMENSION(:) ::   gauss
+#else
+   REAL(wp), INTENT( inout ) :: coeff
+   REAL(wp), POINTER ::   gauss
+#endif
    INTEGER, INTENT( in ) ::  nn_type
    REAL(wp), INTENT( in ) :: rn_sd
@@ -634 +646 @@
    INTEGER :: jklev
 
+#if defined key_traldf_c2d || key_traldf_c3d
    CALL wrk_alloc(jpi,jpj,gauss)
 
@@ -658 +671 @@
        gauss = gauss * rn_sd
        coeff = coeff * ( 1._wp + gauss )
+#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d
        WHERE( coeff > 1._wp ) coeff = 1._wp
        WHERE( coeff < 0._wp ) coeff = 0._wp
+#else
+       IF( coeff > 1._wp ) coeff = 1._wp
+       IF( coeff < 0._wp ) coeff = 0._wp
+#endif
    ELSEIF ( nn_type  == 5 ) THEN
        zsd = rn_sd
@@ -665 +683 @@
        gauss = gauss * zsd + xme
        coeff = exp(gauss) * coeff
+#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d
        WHERE( coeff > 1._wp ) coeff = 1._wp
        WHERE( coeff < 0._wp ) coeff = 0._wp
+#else
+       IF( coeff > 1._wp ) coeff = 1._wp
+       IF( coeff < 0._wp ) coeff = 0._wp
+#endif
    ELSEIF ( nn_type == 6 ) THEN
        zsd = rn_sd
@@ -672 +695 @@
        gauss = gauss * zsd + xme
        coeff = exp(gauss) * coeff
+#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d
        WHERE( coeff > 1._wp ) coeff = 1._wp
        WHERE( coeff < 0._wp ) coeff = 0._wp
+#else
+       IF( coeff > 1._wp ) coeff = 1._wp
+       IF( coeff < 0._wp ) coeff = 0._wp
+#endif
    ELSE
        CALL ctl_stop( 'spp dqdt wrong option')
@@ -687 +715 @@
        jklev = klev
      ELSE
-       jklev = 0 
+       jklev = 0
      ENDIF
      CALL spp_stats(kt,kspp,jklev,coeff)
@@ -694 +722 @@
    CALL wrk_dealloc(jpi,jpj,gauss)
 
+#else
+   CALL ctl_stop( 'spp_gen: parameter perturbation will only work with '// &
+                  'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
    END SUBROUTINE
 
@@ -704 +738 @@
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: mt,kp,kl
-   REAL(wp), INTENT(IN) :: rcf(jpi,jpj)
+#if defined key_traldf_c3d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: rcf
+#elif defined key_traldf_c2d
+   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: rcf
+#elif defined key_traldf_c1d
+   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: rcf
+#else
+   REAL(wp), INTENT( inout ) :: rcf
+#endif
    REAL(wp) :: mi,ma
    CHARACTER(LEN=16) :: cstr = '                '
-   SELECT CASE ( kp ) 
-     CASE( jk_spp_alb   ) 
+   SELECT CASE ( kp )
+     CASE( jk_spp_alb   )
        cstr = 'ALBEDO      '
-     CASE( jk_spp_rhg   ) 
-       cstr = 'ICE RHEOLOGY' 
-     CASE( jk_spp_relw  ) 
-       cstr = 'RELATIVE WND' 
-     CASE( jk_spp_dqdt  ) 
-       cstr = 'SST RELAXAT.' 
-     CASE( jk_spp_deds  ) 
-       cstr = 'SSS RELAXAT.' 
-     CASE( jk_spp_arnf  ) 
-       cstr = 'RIVER MIXING' 
-     CASE( jk_spp_geot  ) 
-       cstr = 'GEOTHERM.FLX' 
-     CASE( jk_spp_qsi0  ) 
-       cstr = 'SOLAR EXTIN.' 
-     CASE( jk_spp_bfr   ) 
-       cstr = 'BOTTOM FRICT' 
-     CASE( jk_spp_aevd  ) 
-       cstr = 'EDDY VISCDIF' 
-     CASE( jk_spp_avt   ) 
+     CASE( jk_spp_rhg   )
+       cstr = 'ICE RHEOLOGY'
+     CASE( jk_spp_relw  )
+       cstr = 'RELATIVE WND'
+     CASE( jk_spp_dqdt  )
+       cstr = 'SST RELAXAT.'
+     CASE( jk_spp_deds  )
+       cstr = 'SSS RELAXAT.'
+     CASE( jk_spp_arnf  )
+       cstr = 'RIVER MIXING'
+     CASE( jk_spp_geot  )
+       cstr = 'GEOTHERM.FLX'
+     CASE( jk_spp_qsi0  )
+       cstr = 'SOLAR EXTIN.'
+     CASE( jk_spp_bfr   )
+       cstr = 'BOTTOM FRICT'
+     CASE( jk_spp_aevd  )
+       cstr = 'EDDY VISCDIF'
+     CASE( jk_spp_avt   )
        cstr = 'VERT. DIFFUS'
-     CASE( jk_spp_avm   ) 
+     CASE( jk_spp_avm   )
        cstr = 'VERT. VISCOS'
-     CASE( jk_spp_tkelc ) 
+     CASE( jk_spp_tkelc )
        cstr = 'TKE LANGMUIR'
-     CASE( jk_spp_tkedf ) 
-       cstr = 'TKE RN_EDIFF' 
-     CASE( jk_spp_tkeds ) 
-       cstr = 'TKE RN_EDISS' 
-     CASE( jk_spp_tkebb ) 
+     CASE( jk_spp_tkedf )
+       cstr = 'TKE RN_EDIFF'
+     CASE( jk_spp_tkeds )
+       cstr = 'TKE RN_EDISS'
+     CASE( jk_spp_tkebb )
        cstr = 'TKE RN_EBB  '
-     CASE( jk_spp_tkefr ) 
+     CASE( jk_spp_tkefr )
        cstr = 'TKE RN_EFR  '
-     CASE( jk_spp_ahtu  ) 
+     CASE( jk_spp_ahtu  )
        cstr = 'TRALDF AHTU '
-     CASE( jk_spp_ahtv  ) 
+     CASE( jk_spp_ahtv  )
        cstr = 'TRALDF AHTV '
-     CASE( jk_spp_ahtw  ) 
+     CASE( jk_spp_ahtw  )
        cstr = 'TRALDF AHTW '
-     CASE( jk_spp_ahtt  ) 
+     CASE( jk_spp_ahtt  )
        cstr = 'TRALDF AHTT '
      CASE( jk_spp_ahubbl )
@@ -765 +807 @@
        CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
    END SELECT
+#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d
    mi = MINVAL(rcf)
    ma = MAXVAL(rcf)
+#else
+   mi = rcf
+   ma = rcf
+#endif
    IF(lk_mpp) CALL mpp_min(mi)
    IF(lk_mpp) CALL mpp_max(ma)
@@ -795 +842 @@
 
    DO jp=1,jk_spp
-     SELECT CASE ( jp ) 
-       CASE( jk_spp_alb   ) 
+     SELECT CASE ( jp )
+       CASE( jk_spp_alb   )
          cstr = 'ALBEDO      '
-       CASE( jk_spp_rhg   ) 
-         cstr = 'ICE RHEOLOGY' 
-       CASE( jk_spp_relw  ) 
-         cstr = 'RELATIVE WND' 
-       CASE( jk_spp_dqdt  ) 
-         cstr = 'SST RELAXAT.' 
-       CASE( jk_spp_deds  ) 
-         cstr = 'SSS RELAXAT.' 
-       CASE( jk_spp_arnf  ) 
-         cstr = 'RIVER MIXING' 
-       CASE( jk_spp_geot  ) 
-         cstr = 'GEOTHERM.FLX' 
-       CASE( jk_spp_qsi0  ) 
-         cstr = 'SOLAR EXTIN.' 
-       CASE( jk_spp_bfr   ) 
-         cstr = 'BOTTOM FRICT' 
-       CASE( jk_spp_aevd  ) 
-         cstr = 'EDDY VISCDIF' 
-       CASE( jk_spp_avt   ) 
+       CASE( jk_spp_rhg   )
+         cstr = 'ICE RHEOLOGY'
+       CASE( jk_spp_relw  )
+         cstr = 'RELATIVE WND'
+       CASE( jk_spp_dqdt  )
+         cstr = 'SST RELAXAT.'
+       CASE( jk_spp_deds  )
+         cstr = 'SSS RELAXAT.'
+       CASE( jk_spp_arnf  )
+         cstr = 'RIVER MIXING'
+       CASE( jk_spp_geot  )
+         cstr = 'GEOTHERM.FLX'
+       CASE( jk_spp_qsi0  )
+         cstr = 'SOLAR EXTIN.'
+       CASE( jk_spp_bfr   )
+         cstr = 'BOTTOM FRICT'
+       CASE( jk_spp_aevd  )
+         cstr = 'EDDY VISCDIF'
+       CASE( jk_spp_avt   )
          cstr = 'VERT. DIFFUS'
-       CASE( jk_spp_avm   ) 
+       CASE( jk_spp_avm   )
          cstr = 'VERT. VISCOS'
-       CASE( jk_spp_tkelc ) 
+       CASE( jk_spp_tkelc )
          cstr = 'TKE LANGMUIR'
-       CASE( jk_spp_tkedf ) 
-         cstr = 'TKE RN_EDIFF' 
-       CASE( jk_spp_tkeds ) 
-         cstr = 'TKE RN_EDISS' 
-       CASE( jk_spp_tkebb ) 
+       CASE( jk_spp_tkedf )
+         cstr = 'TKE RN_EDIFF'
+       CASE( jk_spp_tkeds )
+         cstr = 'TKE RN_EDISS'
+       CASE( jk_spp_tkebb )
          cstr = 'TKE RN_EBB  '
-       CASE( jk_spp_tkefr ) 
+       CASE( jk_spp_tkefr )
          cstr = 'TKE RN_EFR  '
-       CASE( jk_spp_ahtu  ) 
+       CASE( jk_spp_ahtu  )
          cstr = 'TRALDF AHTU '
-       CASE( jk_spp_ahtv  ) 
+       CASE( jk_spp_ahtv  )
          cstr = 'TRALDF AHTV '
-       CASE( jk_spp_ahtw  ) 
+       CASE( jk_spp_ahtw  )
          cstr = 'TRALDF AHTW '
-       CASE( jk_spp_ahtt  ) 
+       CASE( jk_spp_ahtt  )
          cstr = 'TRALDF AHTT '
        CASE( jk_spp_ahubbl )
@@ -1195 +1242 @@
       ! Note: sshn should be staggered before being used.
       SELECT CASE ( cd_type )
-               CASE ( 'T' )  
+               CASE ( 'T' )
                 jk=1
                 zv = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:)*zts(:,:,1) )
@@ -1285 +1332 @@
       ! Random noise on 2d field
       IF ( istep == 1 ) THEN
-         CALL sppt_rand2d( g2d ) 
+         CALL sppt_rand2d( g2d )
          CALL lbc_lnk( g2d , 'T', 1._wp)
          g2d_save = g2d
@@ -1297 +1344 @@
          g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev
       ENDIF
-   
+
       ! Laplacian filter and re-normalization
       DO jf = 1, nk
@@ -1314 +1361 @@
       ENDIF
 #endif
-   
+
       ! AR(1) process and array swap
       g2d = a*gb + b*g2d
@@ -1360 +1407 @@
         ENDDO
       ENDIF
-      
+
       ! Bound
       IF( nn_sppt_step_bound .EQ. 2 ) THEN
@@ -1482 +1529 @@
 
 #ifdef NEMO_V34
-      REWIND( numnam )            
+      REWIND( numnam )
       READ  ( numnam, namstopack )
 #else
-      REWIND( numnam_ref ) 
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp )
 
-      REWIND( numnam_cfg )  
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp )
@@ -1568 +1615 @@
          WRITE(numout,*)
          WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     spp_filter_pass:', spp_filter_pass
-         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev      
+         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev
          WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_spp_tau   :', rn_spp_tau
          WRITE(numout,*)
-         WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr  
-         WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd   
-         WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt 
-         WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd   
-         WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt 
-         WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd   
-         WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt  
-         WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd   
-         WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm  
-         WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd   
+         WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr
+         WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd
+         WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt
+         WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd
+         WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt
+         WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd
+         WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt
+         WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd
+         WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm
+         WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd
          WRITE(numout,*) '       SPP for solar penetration scheme  (only RGB)        nn_spp_qsi0  :', nn_spp_qsi0
-         WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd   
+         WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd
          WRITE(numout,*) '       SPP for horiz. diffusivity  U                       nn_spp_ahtu  :', nn_spp_ahtu
-         WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd
          WRITE(numout,*) '       SPP for horiz. diffusivity  V                       nn_spp_ahtv  :', nn_spp_ahtv
-         WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd
          WRITE(numout,*) '       SPP for horiz. diffusivity  W                       nn_spp_ahtw  :', nn_spp_ahtw
-         WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd
          WRITE(numout,*) '       SPP for horiz. diffusivity  T                       nn_spp_ahtt  :', nn_spp_ahtt
-         WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd
          WRITE(numout,*) '       SPP for horiz. viscosity (1/3)                      nn_spp_ahm1  :', nn_spp_ahm1
-         WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd
          WRITE(numout,*) '       SPP for horiz. viscosity (2/4)                      nn_spp_ahm2  :', nn_spp_ahm2
-         WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd   
+         WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd
          WRITE(numout,*) '       SPP for relative wind factor                        nn_spp_relw  :', nn_spp_relw
          WRITE(numout,*) '       (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]'
-         WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd   
+         WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd
          WRITE(numout,*) '       SPP for mixing close to river mouth                 nn_spp_arnf  :', nn_spp_arnf
-         WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd   
+         WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd
          WRITE(numout,*) '       SPP for geothermal heating                          nn_spp_geot  :', nn_spp_geot
-         WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd   
+         WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd
          WRITE(numout,*) '       SPP for enhanced vertical diffusion                 nn_spp_aevd  :', nn_spp_aevd
-         WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd   
+         WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd
          WRITE(numout,*) '       SPP for TKE rn_lc    Langmuir cell coefficient      nn_spp_tkelc :', nn_spp_tkelc
-         WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd   
+         WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd
          WRITE(numout,*) '       SPP for TKE rn_ediff Eddy diff. coefficient         nn_spp_tkedf :', nn_spp_tkedf
-         WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd   
+         WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd
          WRITE(numout,*) '       SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds
-         WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd   
+         WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd
          WRITE(numout,*) '       SPP for TKE rn_ebb   Surface input of tke           nn_spp_tkebb :', nn_spp_tkebb
-         WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd   
+         WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd
          WRITE(numout,*) '       SPP for TKE rn_efr   Fraction of srf TKE below ML   nn_spp_tkefr :', nn_spp_tkefr
-         WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd   
+         WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd
          WRITE(numout,*) '       SPP for BBL U  diffusivity                          nn_spp_ahubbl:', nn_spp_ahubbl
          WRITE(numout,*) '                                            STDEV          rn_ahubbl_sd :', rn_ahubbl_sd
@@ -1626 +1673 @@
          WRITE(numout,*)
          WRITE(numout,*) ' SKEB Perturbation scheme '
-         WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb    
-         WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb    
+         WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb
+         WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb
          WRITE(numout,*) '       Frequency update for dissipation mask               nn_dcom_freq :', nn_dcom_freq
          WRITE(numout,*) '       Numerical dissipation factor (resolut. dependent)   rn_kh        :', rn_kh
          WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     skeb_filter_pass:', skeb_filter_pass
-         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev      
+         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev
          WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_skeb_tau  :', rn_skeb_tau
          WRITE(numout,*) '       Option of convection energy dissipation             nn_dconv     :', nn_dconv
@@ -1752 +1799 @@
 
       ! Find filter attenuation factor
-   
+
       flt_fac = sppt_fltfac( sppt_filter_pass )
       rdt_sppt = nn_rndm_freq * rn_rdt
-   
+
       IF( ln_sppt_taumap ) THEN
          CALL iom_open ( 'sppt_tau_map', inum )
@@ -1798 +1845 @@
       gauss_b = 0._wp
       ! Weigths
-      gauss_w(:)    = 1.0_wp 
+      gauss_w(:)    = 1.0_wp
       IF( nn_vwei .eq. 1 ) THEN
         gauss_w(1)    = 0.0_wp
@@ -1861 +1908 @@
       IF(lwp .and. ln_stopack_diags) &
       CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-   
+
    END SUBROUTINE stopack_init
    !
@@ -1874 +1921 @@
       INTEGER :: id1, jseed
       CHARACTER(LEN=10)   ::   clseed='spsd0_0000'
-      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type 
-      INTEGER(KIND=8)     ::   ivals(8) 
+      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
+      INTEGER(KIND=8)     ::   ivals(8)
       REAL(wp)            ::   zrseed4(4)           ! RNG seeds in integer type
       REAL(wp)            ::   zrseed2d(jpi,jpj)
@@ -1983 +2030 @@
       !!---------------------------------------------------------------------
       !
-      ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 
-      gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 
+      ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,&
+      gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , &
       spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,&
       gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),&
@@ -2208 +2255 @@
       IF ( ln_dpsiv ) THEN
        DO jp=1,jpni-1
-         IF( jpri == jp ) THEN ! SEND TO EAST 
+         IF( jpri == jp ) THEN ! SEND TO EAST
           zwrk(1:jpj) = dpsiv(jpi-1,:)
           tag=2000+narea
@@ -2268 +2315 @@
    REAL :: ds,dt,dtot,kh2
    INTEGER :: ji,jj,jk
-   
+
    IF ( mt .eq. nit000 ) THEN
         ALLOCATE ( dnum(jpi,jpj,jpk) )
@@ -2287 +2334 @@
           dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e2v(ji,jj) + &
                (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e1u(ji,jj)
-  
+
           dtot = sqrt( ds*ds + dt*dt ) * tmask(ji,jj,jk)
           dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(ji,jj)*e2t(ji,jj)
@@ -2293 +2340 @@
       ENDDO
      ENDDO
-   
+
      CALL lbc_lnk(dnum,'T',1._wp)
 
    ENDIF
 
-   END SUBROUTINE 
+   END SUBROUTINE
 
    SUBROUTINE skeb_dcon ( mt )
@@ -2329 +2376 @@
 
            zz = - grav*avt(ji,jj,jk) * ( rhd(ji,jj,jk)-rhd(ji,jj,jk-1) ) * wmask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) &
-              & / ( rau0 * fse3w(ji,jj,jk) ) 
+              & / ( rau0 * fse3w(ji,jj,jk) )
 
            dcon(ji,jj,jk) = kc2*zz*e1t(ji,jj)*e2t(ji,jj)*rau0 / fse3w(ji,jj,jk)
@@ -2378 +2425 @@
      IF(ln_skeb_own_gauss) THEN
        DO jk=1,jpkm1
-         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 
+         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:)
        ENDDO
      ELSE
@@ -2407 +2454 @@
      IF(ln_skeb_own_gauss) THEN
        DO jk=1,jpkm1
-         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
+         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
        ENDDO
      ELSE
@@ -2440 +2487 @@
    IF(ln_skeb_own_gauss) THEN
      DO jk=1,jpkm1
-       psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
+       psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
      ENDDO
    ELSE

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

@@ -12 +12 @@
 
    !!----------------------------------------------------------------------
-   !!   tra_bbc      : update the tracer trend at ocean bottom 
+   !!   tra_bbc      : update the tracer trend at ocean bottom
    !!   tra_bbc_init : initialization of geothermal heat flux trend
    !!----------------------------------------------------------------------
@@ -19 +19 @@
    USE phycst          ! physical constants
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager
@@ -44 +44 @@
    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read)
- 
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -58 +58 @@
       !!                  ***  ROUTINE tra_bbc  ***
       !!
-      !! ** Purpose :   Compute the bottom boundary contition on temperature 
-      !!              associated with geothermal heating and add it to the 
+      !! ** Purpose :   Compute the bottom boundary contition on temperature
+      !!              associated with geothermal heating and add it to the
       !!              general trend of temperature equations.
       !!
-      !! ** Method  :   The geothermal heat flux set to its constant value of 
+      !! ** Method  :   The geothermal heat flux set to its constant value of
       !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999).
       !!       The temperature trend associated to this heat flux through the
@@ -92 +92 @@
       !                             !  Add the geothermal heat flux trend on temperature
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_geot > 0) THEN
           qgh_trd1(:,:) = qgh_trd0(:,:)
           CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot)
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_geot > 0 ) &
+         & CALL ctl_stop( 'tra_bbc: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
@@ -144 +152 @@
       CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files
       !
-      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
+      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir
       !!----------------------------------------------------------------------
 

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

@@ -32 +32 @@
    USE trdtra         ! trends: active tracers
    !
-   USE iom            ! IOM library               
+   USE iom            ! IOM library
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions
@@ -38 +38 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE stopack
 
@@ -199 +199 @@
       !
       ahu_bbl_1(:,:) = ahu_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
           CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
+
       ahv_bbl_1(:,:) = ahv_bbl(:,:)
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
           CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
       ENDIF
+#else
+      IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) &
+         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
+
       !
       DO jn = 1, kjpt                                     ! tracer loop
@@ -215 +230 @@
             END DO
          END DO
-         !               
+         !
          DO jj = 2, jpjm1                                    ! Compute the trend
             DO ji = 2, jpim1
@@ -426 +441 @@
                   za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
                   zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
-                  !                                                          ! 2*masked bottom density gradient 
+                  !                                                          ! 2*masked bottom density gradient
                   zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
                             - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
@@ -585 +600 @@
                gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
             ENDIF
-            !     
+            !
             IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
                mgrhv(ji,jj) = INT(  SIGN( 1.e0, &

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -9 +9 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
-   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
-   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
+   !!            3.2  !  2009-04  (G. Madec & NEMO team)
+   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model
    !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
    !!----------------------------------------------------------------------
@@ -43 +43 @@
    !                                 !!* Namelist namtra_qsr: penetrative solar radiation
    LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag
-   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag  
+   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
@@ -51 +51 @@
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
- 
+
    ! Module variables
    REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0
@@ -80 +80 @@
       !!      Considering the 2 wavebands case:
       !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
-      !!         The temperature trend associated with the solar radiation penetration 
+      !!         The temperature trend associated with the solar radiation penetration
       !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
       !!         At the bottom, boudary condition for the radiation is no flux :
@@ -86 +86 @@
       !!      in the last ocean level.
       !!         In z-coordinate case, the computation is only done down to the
-      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients 
+      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
       !!      used for the computation are calculated one for once as they
       !!      depends on k only.
@@ -106 +106 @@
       REAL(wp) ::   zz0, zz1, z1_e3t     !    -         -
       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
-      REAL(wp) ::   zlogc, zlogc2, zlogc3 
+      REAL(wp) ::   zlogc, zlogc2, zlogc3
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d
@@ -113 +113 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( kt == nit000 ) THEN
@@ -124 +124 @@
 
       IF( l_trdtra ) THEN      ! Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
       ENDIF
@@ -150 +150 @@
       !                                        Compute now qsr tracer content field
       !                                        ************************************
-      
+
       !                                           ! ============================================== !
       IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
@@ -159 +159 @@
          !                                        Add to the general trend
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -180 +180 @@
          ENDIF
          !                                        ! ============================================== !
-      ELSE                                        !  Ocean alone : 
+      ELSE                                        !  Ocean alone :
          !                                        ! ============================================== !
          !
-         ! 
+         !
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN
-             xsi0r = rn_si0
-             CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 )
-             xsi0r = 1.e0 / xsi0r
-         ENDIF
+            xsi0r = rn_si0
+            CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 )
+            xsi0r = 1.e0 / xsi0r
+         ENDIF
+#else
+         IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) &
+            & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
+
          !                                               ! ------------------------- !
          IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
@@ -199 +206 @@
                   CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
+                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)
                   ENDDO
                   !
@@ -220 +227 @@
                         zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
                         zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
-                        zCze    = 1.12  * (zchl)**0.803 
+                        zCze    = 1.12  * (zchl)**0.803
                         DO jk = 1, nksr + 1
                            zpsi = fsdept(ji,jj,jk) / zze
@@ -231 +238 @@
                ELSE                              !* Variable ocean volume but constant chrlorophyll
                   DO jk = 1, nksr + 1
-                     zchl3d(:,:,jk) = 0.05 
+                     zchl3d(:,:,jk) = 0.05
                   ENDDO
                ENDIF
@@ -256 +263 @@
 !CDIR NOVERRCHK
                   DO jj = 1, jpj
-!CDIR NOVERRCHK   
+!CDIR NOVERRCHK
                      DO ji = 1, jpi
                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) )
@@ -289 +296 @@
                      END DO
                   END DO
-                  ! 
+                  !
                   DO jj = 1, jpj
                      DO ji = 1, jpi
@@ -296 +303 @@
                         zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
                         zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
-                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
+                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)
                      END DO
                   END DO
@@ -321 +328 @@
                zz0   =        rn_abs   * r1_rau0_rcp
                zz1   = ( 1. - rn_abs ) * r1_rau0_rcp
-               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
+               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
                   DO jj = 1, jpj
                      DO ji = 1, jpi
                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
+                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )
                      END DO
                   END DO
@@ -344 +351 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        ! (ISF) no light penetration below the ice shelves         
+                        ! (ISF) no light penetration below the ice shelves
                         qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1)
                      END DO
@@ -360 +367 @@
          !                                        Add to the general trend
          DO jk = 1, nksr
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   z1_e3t = zfact / fse3t(ji,jj,jk)
@@ -376 +383 @@
             &                    'at it= ', kt,' date= ', ndastp
          IF(lwp) WRITE(numout,*) '~~~~'
-         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
+         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
-         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
+         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm
          IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
          !
@@ -386 +393 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )
       ENDIF
       !                       ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr')
@@ -408 +415 @@
       !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
       !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
-      !!      default values correspond to clear water (type I in Jerlov' 
+      !!      default values correspond to clear water (type I in Jerlov'
       !!      (1968) classification.
       !!         called by tra_qsr at the first timestep (nit000)
@@ -435 +442 @@
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
       !
-      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
 
@@ -465 +472 @@
 
       IF( ln_traqsr ) THEN     ! control consistency
-         !                      
+         !
          IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN
             CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
@@ -480 +487 @@
             &              ' 2 bands, 3 RGB bands or bio-model light penetration' )
          !
-         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
+         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1
          IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2
          IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3
@@ -497 +504 @@
       ENDIF
       !                          ! ===================================== !
-      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !  
+      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !
          !                       ! ===================================== !
          !
@@ -539 +546 @@
                   zchl = 0.05                                 ! constant chlorophyll
                   irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll 
+                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll
                   zekg(:,:) = rkrgb(2,irgb)
                   zekr(:,:) = rkrgb(3,irgb)
@@ -546 +553 @@
                   ze0(:,:,1) = rn_abs
                   ze1(:,:,1) = zcoef
-                  ze2(:,:,1) = zcoef 
+                  ze2(:,:,1) = zcoef
                   ze3(:,:,1) = zcoef
                   zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask
-               
+
                   DO jk = 2, nksr+1
 !CDIR NOVERRCHK
                      DO jj = 1, jpj
-!CDIR NOVERRCHK   
+!CDIR NOVERRCHK
                         DO ji = 1, jpi
                            zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) )
@@ -566 +573 @@
                         END DO
                      END DO
-                  END DO 
+                  END DO
                   !
                   DO jk = 1, nksr
@@ -598 +605 @@
                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
-                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
+                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)
                      END DO
                   END DO
@@ -607 +614 @@
          ENDIF
          !                       ! ===================================== !
-      ELSE                       !        No light penetration           !                   
+      ELSE                       !        No light penetration           !
          !                       ! ===================================== !
          IF(lwp) THEN
@@ -625 +632 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
+      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
+      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init')

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

@@ -107 +107 @@
          WRITE(numout,*) '~~~~~~~~'
       ENDIF
-      ! 
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN
          bfrcoef2d = bfrcoef2d0
          CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr)
       ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) &
+         & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
@@ -140 +146 @@
                END DO
             END IF
-         !   
+         !
          ELSE
             zbfrt(:,:) = bfrcoef2d(:,:)
@@ -183 +189 @@
                DO ji = 2, jpim1
                   ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
                   ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
                   !
@@ -363 +369 @@
             bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
          ENDIF
-         
+
          IF ( ln_isfcav ) THEN
             IF(ln_tfr2d) THEN

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

@@ -20 +20 @@
    USE zdfkpp          ! KPP vertical mixing
    USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   USE trdtra          ! trends manager: tracers
    USE in_out_manager  ! I/O manager
    USE iom             ! for iom_put
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE timing          ! Timing
-   USE stopack         
+   USE stopack
 
    IMPLICIT NONE
@@ -45 +45 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_evd  ***
-      !!                   
+      !!
       !! ** Purpose :   Local increased the vertical eddy viscosity and diffu-
       !!      sivity coefficients when a static instability is encountered.
       !!
       !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients
-      !!      are set to avevd (namelist parameter) if the water column is 
+      !!      are set to avevd (namelist parameter) if the water column is
       !!      statically unstable (i.e. if rn2 < -1.e-12 )
       !!
@@ -77 +77 @@
       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
 
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN
          rn_avevd0(:,:) = rn_avevd
          CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd)
       ENDIF
+#else
+      IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) &
+         & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
 
       SELECT CASE ( nn_evdm )
@@ -88 +94 @@
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
-         DO jk = 1, jpkm1 
+         DO jk = 1, jpkm1
             DO jj = 2, jpj             ! no vector opt.
                DO ji = 2, jpi
@@ -106 +112 @@
                END DO
             END DO
-         END DO 
+         END DO
          CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
@@ -113 +119 @@
          CALL iom_put( "avm_evd", zavm_evd )              ! output this change
          !
-      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12)
          DO jk = 1, jpkm1
-!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
+!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!
             DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
                DO ji = 1, jpi
 #if defined key_zdfkpp
                   ! no evd mixing in the boundary layer with KPP
-                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &          
+                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &
 #else
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
@@ -129 +135 @@
          END DO
          !
-      END SELECT 
+      END SELECT
 
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdfgls  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the gls 
+   !! Ocean physics:  vertical mixing coefficient computed from the gls
    !!                 turbulent closure parameterization
    !!======================================================================
@@ -16 +16 @@
    !!   gls_rst       : read/write gls restart in ocean restart file
    !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers 
+   USE oce            ! ocean dynamics and active tracers
    USE dom_oce        ! ocean space and time domain
    USE domvvl         ! ocean space and time domain : variable volume layer
@@ -31 +31 @@
    USE iom            ! I/O manager library
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE stopack
 
@@ -62 +62 @@
    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
-   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
-   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn    
+   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1
+   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn
    REAL(wp) ::   rghmin        = -0.28_wp
    REAL(wp) ::   rgh0          =  0.0329_wp
@@ -73 +73 @@
    REAL(wp) ::   ra2           =  0.74_wp
    REAL(wp) ::   rb1           = 16.60_wp
-   REAL(wp) ::   rb2           = 10.10_wp         
-   REAL(wp) ::   re2           =  1.33_wp         
+   REAL(wp) ::   rb2           = 10.10_wp
+   REAL(wp) ::   re2           =  1.33_wp
    REAL(wp) ::   rl1           =  0.107_wp
    REAL(wp) ::   rl2           =  0.0032_wp
@@ -134 +134 @@
       INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
       REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
-      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
+      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -
       REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
@@ -140 +140 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
@@ -157 +157 @@
       CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )
       CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
-      
+
       ! Preliminary computing
 
@@ -166 +166 @@
          avm (:,:,:) = avm_k (:,:,:)
          avmu(:,:,:) = avmu_k(:,:,:)
-         avmv(:,:,:) = avmv_k(:,:,:) 
+         avmv(:,:,:) = avmv_k(:,:,:)
       ENDIF
 
       ! Compute surface and bottom friction at T-points
-!CDIR NOVERRCHK          
-      DO jj = 2, jpjm1          
-!CDIR NOVERRCHK         
-         DO ji = fs_2, fs_jpim1   ! vector opt.         
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
             !
             ! surface friction
             ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
-            !   
-            ! bottom friction (explicit before friction)        
-            ! Note that we chose here not to bound the friction as in dynbfr)   
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
-            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
-         END DO         
-      END DO    
+            !
+            ! bottom friction (explicit before friction)
+            ! Note that we chose here not to bound the friction as in dynbfr)
+            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &
+               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )
+            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &
+               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )
+            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+         END DO
+      END DO
 
       ! Set surface roughness length
       SELECT CASE ( nn_z0_met )
       !
-      CASE ( 0 )             ! Constant roughness          
+      CASE ( 0 )             ! Constant roughness
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
@@ -227 +227 @@
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1 
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)
@@ -276 +276 @@
                IF( ln_sigpsi ) THEN
                   zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
-                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
+                  zwall_psi(ji,jj,jk) = rsc_psi /   &
                      &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
                ELSE
@@ -295 +295 @@
                ! diagonal
                z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
+                  &                       + rdt * zdiss * tmask(ji,jj,jk)
                !
                ! right hand side in en
@@ -317 +317 @@
       ! First level
       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
+      en(:,:,1) = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
       z_elem_b(:,:,1) = 1._wp
-      ! 
+      !
       ! One level below
       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
+      z_elem_a(:,:,2) = 0._wp
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
@@ -335 +335 @@
       ! Dirichlet conditions at k=1
       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
+      en(:,:,1)       = MAX(en(:,:,1), rn_emin)
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0._wp
@@ -358 +358 @@
       SELECT CASE ( nn_bc_bot )
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
@@ -378 +378 @@
                z_elem_c(ji,jj,ibotm1) = 0._wp
                z_elem_b(ji,jj,ibotm1) = 1._wp
-               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
+               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
             END DO
          END DO
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -429 +429 @@
          END DO
       END DO
-      !                                            ! set the minimum value of tke 
+      !                                            ! set the minimum value of tke
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
 
@@ -471 +471 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
+                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn
                END DO
             END DO
@@ -490 +490 @@
                !
                ! psi / k
-               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
+               zratio = psi(ji,jj,jk) / eb(ji,jj,jk)
                !
                ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
@@ -510 +510 @@
                zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
                zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
-               !                                                        
+               !
                ! building the matrix
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
@@ -552 +552 @@
       z_elem_c(:,:,2) = 0._wp
       z_elem_b(:,:,2) = 1._wp
-      ! 
+      !
       !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
@@ -576 +576 @@
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
 
-      !   
+      !
       !
       END SELECT
@@ -586 +586 @@
       !
       !
-      CASE ( 0 )             ! Dirichlet 
+      CASE ( 0 )             ! Dirichlet
          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
@@ -611 +611 @@
          !
       CASE ( 1 )             ! Neumman boundary condition
-         !                      
+         !
 !CDIR NOVERRCHK
          DO jj = 2, jpjm1
@@ -693 +693 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
+                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)
                END DO
             END DO
@@ -720 +720 @@
                eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
                mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
-               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
+               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
             END DO
          END DO
-      END DO 
+      END DO
 
       !
@@ -807 +807 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack) THEN
             IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk)
             IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk)
-         ENDIF 
+         ENDIF
+#else
+         IF ( ln_stopack ) &
+            & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       !
@@ -851 +857 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_gls_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!      viscosity when using a gls turbulent closure scheme
       !!
@@ -1071 +1077 @@
          !
       END SELECT
-    
+
       !                                !* Set Schmidt number for psi diffusion in the wave breaking case
       !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009
       !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001
       IF( ln_sigpsi ) THEN
-         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
+         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf
          ! Verification: retrieve Burchard (2001) results by uncomenting the line below:
          ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work
@@ -1084 +1090 @@
          rsc_psi0 = rsc_psi
       ENDIF
- 
+
       !                                !* Shear free turbulence parameters
       !
@@ -1130 +1136 @@
       rc04  = rc03 * rc0
       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
+      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking
       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
-      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
+      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness
       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
-      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
+      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking
 
       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
@@ -1151 +1157 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      !                              
+      !
       CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
       !
@@ -1162 +1168 @@
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE ts_rst  ***
-      !!                     
+      !!
       !! ** Purpose :   Read or write TKE file (en) in restart file
       !!
       !! ** Method  :   use of IOM library
-      !!                if the restart does not contain TKE, en is either 
+      !!                if the restart does not contain TKE, en is either
       !!                set to rn_emin or recomputed (nn_igls/=0)
       !!----------------------------------------------------------------------
@@ -1178 +1184 @@
       !!----------------------------------------------------------------------
       !
-      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
          !                                   ! ---------------
          IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1197 +1203 @@
                CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
                IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
-            ELSE                        
+            ELSE
                IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
                en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.05        
+               mxln(:,:,:) = 0.05
                avt_k (:,:,:) = avt (:,:,:)
                avm_k (:,:,:) = avm (:,:,:)
@@ -1210 +1216 @@
             IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            mxln(:,:,:) = 0.05
          ENDIF
          !
@@ -1217 +1223 @@
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
          IF(nn_timing == 2)  CALL timing_start('iom_rstput')
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
          CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
          CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
          CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
          CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
@@ -1239 +1245 @@
       WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
    END SUBROUTINE zdf_gls_init
-   
+
    SUBROUTINE zdf_gls( kt )          ! Empty routine
    IMPLICIT NONE
-      INTEGER, INTENT(in) ::   kt ! ocean time step   
+      INTEGER, INTENT(in) ::   kt ! ocean time step
       WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_gls
-   
+
    SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
    IMPLICIT NONE
@@ -1256 +1262 @@
    !!======================================================================
 END MODULE zdfgls
-

branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -2 +2 @@
    !!======================================================================
    !!                       ***  MODULE  zdftke  ***
-   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !! Ocean physics:  vertical mixing coefficient computed from the tke
    !!                 turbulent closure parameterization
    !!=====================================================================
@@ -22 +22 @@
    !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
    !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
-   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
-   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
@@ -52 +52 @@
    USE wrk_nemo       ! work arrays
    USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
 #if defined key_agrif
    USE agrif_opa_interp
@@ -74 +74 @@
    INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
    REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
-   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
    REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
    REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
@@ -102 +102 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4
-   
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -118 +118 @@
       !!----------------------------------------------------------------------
       ALLOCATE(                                                                    &
-         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &      
+         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &
 #if defined key_c1d
          &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
@@ -147 +147 @@
       !!         surface: en = max( rn_emin0, rn_ebb * taum )
       !!         bottom : en = rn_emin
-      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
-      !!
-      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following
       !!      time stepping: implicit for vertical diffusion term, linearized semi
-      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
-      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal
       !!      linear system is solved. Note that buoyancy and shear terms are
       !!      discretized in a energy conserving form (Bruchard 2002).
@@ -160 +160 @@
       !!
       !!        The now vertical eddy vicosity and diffusivity coefficients are
-      !!      given by: 
+      !!      given by:
       !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!              avt = max( avmb, pdl * avm                 )  
+      !!              avt = max( avmb, pdl * avm                 )
       !!              eav = max( avmb, avm )
       !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
-      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
       !!
       !! ** Action  :   compute en (now turbulent kinetic energy)
@@ -180 +180 @@
       !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:) 
-         avm (:,:,:) = avm_k (:,:,:) 
-         avmu(:,:,:) = avmu_k(:,:,:) 
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF
-      !
+         avt (:,:,:) = avt_k (:,:,:)
+         avm (:,:,:) = avm_k (:,:,:)
+         avmu(:,:,:) = avmu_k(:,:,:)
+         avmv(:,:,:) = avmv_k(:,:,:)
+      ENDIF
+      !
+#if defined key_traldf_c2d || key_traldf_c3d
       IF( ln_stopack) THEN
          IF( nn_spp_tkelc > 0 ) THEN
@@ -208 +209 @@
          ENDIF
       ENDIF
+#else
+      IF ( ln_stopack ) &
+         & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// &
+                          'key_traldf_c2d or key_traldf_c3d')
+#endif
       !
       CALL tke_tke      ! now tke (en)
@@ -213 +219 @@
       CALL tke_avn      ! now avt, avm, avmu, avmv
       !
-      avt_k (:,:,:) = avt (:,:,:) 
-      avm_k (:,:,:) = avm (:,:,:) 
-      avmu_k(:,:,:) = avmu(:,:,:) 
-      avmv_k(:,:,:) = avmv(:,:,:) 
+      avt_k (:,:,:) = avt (:,:,:)
+      avm_k (:,:,:) = avm (:,:,:)
+      avmu_k(:,:,:) = avmu(:,:,:)
+      avmv_k(:,:,:) = avmv(:,:,:)
       !
 #if defined key_agrif
-      ! Update child grid f => parent grid 
+      ! Update child grid f => parent grid
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
-#endif      
+#endif
       IF (  kt == nitend ) THEN
         DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 )
@@ -271 +277 @@
       CALL wrk_alloc( jpi,jpj, zfact2 )
       CALL wrk_alloc( jpi,jpj, zfact3 )
-      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
       !
       zbbrau = rn_ebb0 / rau0       ! Local constant initialisation
-      zfact1 = -.5_wp * rdt 
+      zfact1 = -.5_wp * rdt
       zfact2 = 1.5_wp * rdt * rn_ediss0
       zfact3 = 0.5_wp       * rn_ediss0
@@ -294 +300 @@
          END DO
       END DO
-      
+
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -333 +339 @@
          imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
          DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
                DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
                   zus  = zcof * taum(ji,jj)
@@ -341 +347 @@
          END DO
          !                               ! finite LC depth
-         DO jj = 1, jpj 
+         DO jj = 1, jpj
             DO ji = 1, jpi
                zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
@@ -378 +384 @@
             DO ji = 1, jpi
                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
                   &                            / (  fse3uw_n(ji,jj,jk)               &
                   &                              *  fse3uw_b(ji,jj,jk)  )
@@ -407 +413 @@
                   !                                                           ! shear prod. at w-point weightened by mask
                zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
                   !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
@@ -464 +470 @@
       END DO
 
-      !                                 ! Save TKE prior to nn_etau addition  
-      e_niw(:,:,:) = en(:,:,:)  
-      !  
+      !                                 ! Save TKE prior to nn_etau addition
+      e_niw(:,:,:) = en(:,:,:)
+      !
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  TKE due to surface and internal wave breaking
@@ -508 +514 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
-                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
+                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
@@ -517 +523 @@
          END DO
       ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up)
-         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
+         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -535 +541 @@
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term  
-         DO jj = 2, jpjm1  
-            DO ji = fs_2, fs_jpim1   ! vector opt.  
-               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)  
-            END DO  
-         END DO  
-      END DO  
-      !  
-      CALL lbc_lnk( e_niw, 'W', 1. )  
+      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( e_niw, 'W', 1. )
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_dealloc( jpi,jpj, zhlc ) 
-      CALL wrk_dealloc( jpi,jpj, zbbrau ) 
-      CALL wrk_dealloc( jpi,jpj, zfact2 ) 
-      CALL wrk_dealloc( jpi,jpj, zfact3 ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_dealloc( jpi,jpj, zhlc )
+      CALL wrk_dealloc( jpi,jpj, zbbrau )
+      CALL wrk_dealloc( jpi,jpj, zfact2 )
+      CALL wrk_dealloc( jpi,jpj, zfact3 )
+      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -563 +569 @@
       !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
       !!
-      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
-      !!              the tke_tke routine). First, the now mixing lenth is 
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
+      !!              the tke_tke routine). First, the now mixing lenth is
       !!      computed from en and the strafification (N^2), then the mixings
       !!      coefficients are computed.
       !!              - Mixing length : a first evaluation of the mixing lengh
       !!      scales is:
-      !!                      mxl = sqrt(2*en) / N  
+      !!                      mxl = sqrt(2*en) / N
       !!      where N is the brunt-vaisala frequency, with a minimum value set
       !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
-      !!        The mixing and dissipative length scale are bound as follow : 
+      !!        The mixing and dissipative length scale are bound as follow :
       !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
       !!                        zmxld = zmxlm = mxl
       !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
-      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
       !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
       !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
-      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
-      !!                    the surface to obtain ldown. the resulting length 
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
+      !!                    the surface to obtain ldown. the resulting length
       !!                    scales are:
-      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxld = sqrt( lup * ldown )
       !!                        zmxlm = min ( lup , ldown )
       !!              - Vertical eddy viscosity and diffusivity:
       !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
-      !!                      avt = max( avmb, pdlr * avm )  
+      !!                      avt = max( avmb, pdlr * avm )
       !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
       !!
@@ -601 +607 @@
       IF( nn_timing == 1 )  CALL timing_start('tke_avn')
 
-      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
 
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -610 +616 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
-      zmxlm(:,:,:)  = rmxl_min    
+      zmxlm(:,:,:)  = rmxl_min
       zmxld(:,:,:)  = rmxl_min
       !
@@ -620 +626 @@
             END DO
          END DO
-      ELSE 
+      ELSE
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
@@ -638 +644 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
+      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
       !
@@ -738 +744 @@
             END DO
          END DO
+#if defined key_traldf_c2d || key_traldf_c3d
          IF( ln_stopack) THEN
             IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk)
             IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk)
          ENDIF
+#else
+         IF ( ln_stopack  ) &
+            & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
       END DO
       CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
@@ -787 +799 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
@@ -797 +809 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
       !!              viscosity when using a tke turbulent closure scheme
       !!
@@ -814 +826 @@
          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
          &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
-         &                 nn_etau , nn_htau  , rn_efr , rn_c   
+         &                 nn_etau , nn_htau  , rn_efr , rn_c
       !!----------------------------------------------------------------------
 
@@ -884 +896 @@
       ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb
       ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr
-      
+
       IF( nn_etau == 2  ) THEN
           ierr = zdf_mxl_alloc()
@@ -896 +908 @@
 
       !                               !* depth of penetration of surface tke
-      IF( nn_etau /= 0 ) THEN      
+      IF( nn_etau /= 0 ) THEN
          htau(:,:) = 0._wp
          SELECT CASE( nn_htau )             ! Choice of the depth of penetration
@@ -902 +914 @@
             htau(:,:) = 10._wp
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
-            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
          CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
             rhtau = -1._wp / LOG( 1._wp - rn_c )
@@ -945 +957 @@
       END DO
       dissl(:,:,:) = 1.e-12_wp
-      !                              
+      !
       CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
       !
@@ -954 +966 @@
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE tke_rst  ***
-     !!                     
+     !!
      !! ** Purpose :   Read or write TKE file (en) in restart file
      !!
      !! ** Method  :   use of IOM library
-     !!                if the restart does not contain TKE, en is either 
-     !!                set to rn_emin or recomputed 
+     !!                if the restart does not contain TKE, en is either
+     !!                set to rn_emin or recomputed
      !!----------------------------------------------------------------------
      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
@@ -968 +980 @@
      !!----------------------------------------------------------------------
      !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
         !                                   ! ---------------
         IF( ln_rstart ) THEN                   !* Read the restart file
@@ -1006 +1018 @@
         avmu_k(:,:,:) = avmu(:,:,:)
         avmv_k(:,:,:) = avmv(:,:,:)
-        
+
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------

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

@@ -110 +110 @@
       IF( lk_tide    )   CALL sbc_tide( kstp )
       IF( lk_bdy     )  THEN
-         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
+         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib
                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
       ENDIF
@@ -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)
@@ -153 +153 @@
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
+#if defined key_traldf_c2d || key_traldf_c3d
          IF ( ln_stopack .AND. ( nn_spp_arnf > 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
+#else
+         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) &
+            & CALL ctl_stop( 'stp: parameter perturbation will only work with '// &
+                             'key_traldf_c2d or key_traldf_c3d')
+#endif
          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
       ENDIF
@@ -204 +210 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
-      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
-                         CALL wzv           ( kstp )  ! now cross-level velocity 
-
-      IF( lk_dynspg_ts ) THEN 
+      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
+                         CALL wzv           ( kstp )  ! now cross-level velocity
+
+      IF( lk_dynspg_ts ) THEN
           ! In case the time splitting case, update almost all momentum trends here:
           ! Note that the computation of vertical velocity above, hence "after" sea level
@@ -244 +250 @@
                                   CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
           IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
-                                  CALL wzv           ( kstp )  ! now cross-level velocity 
+                                  CALL wzv           ( kstp )  ! now cross-level velocity
       ENDIF
 
@@ -320 +326 @@
                &             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)                                   & 
+         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 ,   &
@@ -334 +340 @@
                                 & 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)                                   & 
+                                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 ,   &
@@ -349 +355 @@
                                ua(:,:,:) = ua_sv(:,:,:)
                                va(:,:,:) = va_sv(:,:,:)
-                                                             ! Revert now divergence and rotational to previously computed ones 
+                                                             ! Revert now divergence and rotational to previously computed ones
                                                              !(needed because of the time swap in div_cur, at the beginning of each time step)
                                hdivn(:,:,:) = hdivb(:,:,:)
-                               rotn(:,:,:)  = rotb(:,:,:) 
+                               rotn(:,:,:)  = rotb(:,:,:)
 
                                CALL dyn_bfr( kstp )         ! bottom friction
@@ -398 +404 @@
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! AGRIF
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
-                               CALL Agrif_Integrate_ChildGrids( stp )  
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                               CALL Agrif_Integrate_ChildGrids( stp )
 
       IF ( Agrif_NbStepint().EQ.0 ) THEN
@@ -431 +437 @@
       !
 #if defined key_iomput
-      IF( kstp == nitend .OR. indic < 0 ) THEN 
+      IF( kstp == nitend .OR. indic < 0 ) THEN
                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
-         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
       ENDIF
 #endif
       !
       IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
-      !     
+      !
       !
    END SUBROUTINE stp