Changeset 13355


Ignore:
Timestamp:
2020-07-30T12:12:41+02:00 (12 days ago)
Author:
jwhile
Message:

Merged in changes to fix 1d running - documented in UKMO ocean utils ticket 367

Location:
branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

    r13316 r13355  
    3535   USE oce, ONLY:           & ! active tracer variables 
    3636      & tsn 
    37    USE zdfmxl, ONLY :       & ! mixed layer depth  
     37   USE zdfmxl, ONLY :       & ! mixed layer depth 
    3838#if defined key_karaml 
    3939      & hmld_kara,          & 
    4040      & ln_kara,            & 
    41 #endif    
    42       & hmld,               &  
     41#endif 
     42      & hmld,               & 
    4343      & hmlp,               & 
    44       & hmlpt  
     44      & hmlpt 
    4545   USE asmpar, ONLY:        & ! assimilation parameters 
    4646      & c_asmbkg,           & 
     
    9797 
    9898   IMPLICIT NONE 
    99    PRIVATE                    
     99   PRIVATE 
    100100 
    101101   PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90 
     
    304304 
    305305      ! Allocate and read increments 
    306        
     306 
    307307      IF ( ln_slchltotinc ) THEN 
    308308         ALLOCATE( slchltot_bkginc(jpi,jpj) ) 
    309309         CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc ) 
    310310      ENDIF 
    311        
     311 
    312312      IF ( ln_slchldiainc ) THEN 
    313313         ALLOCATE( slchldia_bkginc(jpi,jpj) ) 
    314314         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc ) 
    315315      ENDIF 
    316        
     316 
    317317      IF ( ln_slchlnoninc ) THEN 
    318318         ALLOCATE( slchlnon_bkginc(jpi,jpj) ) 
    319319         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc ) 
    320320      ENDIF 
    321        
     321 
    322322      IF ( ln_schltotinc ) THEN 
    323323         ALLOCATE( schltot_bkginc(jpi,jpj) ) 
    324324         CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc ) 
    325325      ENDIF 
    326        
     326 
    327327      IF ( ln_slphytotinc ) THEN 
    328328         ALLOCATE( slphytot_bkginc(jpi,jpj) ) 
    329329         CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc ) 
    330330      ENDIF 
    331        
     331 
    332332      IF ( ln_slphydiainc ) THEN 
    333333         ALLOCATE( slphydia_bkginc(jpi,jpj) ) 
    334334         CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc ) 
    335335      ENDIF 
    336        
     336 
    337337      IF ( ln_slphynoninc ) THEN 
    338338         ALLOCATE( slphynon_bkginc(jpi,jpj) ) 
     
    349349         CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc ) 
    350350      ENDIF 
    351        
     351 
    352352      IF ( ln_plchltotinc ) THEN 
    353353         ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) ) 
    354354         CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc ) 
    355355      ENDIF 
    356        
     356 
    357357      IF ( ln_pchltotinc ) THEN 
    358358         ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) ) 
    359359         CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc ) 
    360360      ENDIF 
    361        
     361 
    362362      IF ( ln_pno3inc ) THEN 
    363363         ALLOCATE( pno3_bkginc(jpi,jpj,jpk) ) 
    364364         CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc ) 
    365365      ENDIF 
    366        
     366 
    367367      IF ( ln_psi4inc ) THEN 
    368368         ALLOCATE( psi4_bkginc(jpi,jpj,jpk) ) 
    369369         CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc ) 
    370370      ENDIF 
    371        
     371 
    372372      IF ( ln_pdicinc ) THEN 
    373373         ALLOCATE( pdic_bkginc(jpi,jpj,jpk) ) 
    374374         CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc ) 
    375375      ENDIF 
    376        
     376 
    377377      IF ( ln_palkinc ) THEN 
    378378         ALLOCATE( palk_bkginc(jpi,jpj,jpk) ) 
    379379         CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc ) 
    380380      ENDIF 
    381        
     381 
    382382      IF ( ln_pphinc ) THEN 
    383383         ALLOCATE( pph_bkginc(jpi,jpj,jpk) ) 
    384384         CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc ) 
    385385      ENDIF 
    386        
     386 
    387387      IF ( ln_po2inc ) THEN 
    388388         ALLOCATE( po2_bkginc(jpi,jpj,jpk) ) 
     
    391391 
    392392      ! Allocate balancing increments 
    393        
     393 
    394394      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
    395395         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
     
    402402#endif 
    403403      ENDIF 
    404        
     404 
    405405      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN 
    406406#if defined key_top 
     
    457457      ! Initialise 
    458458      p_incs(:,:) = 0.0 
    459        
     459 
    460460      ! read from file 
    461461      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 ) 
    462        
     462 
    463463      ! Apply the masks 
    464464      p_incs(:,:) = p_incs(:,:) * tmask(:,:,1) 
    465        
     465 
    466466      ! Set missing increments to 0.0 rather than 1e+20 
    467467      ! to allow for differences in masks 
     
    495495      ! Initialise 
    496496      p_incs(:,:,:) = 0.0 
    497        
     497 
    498498      ! read from file 
    499499      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 ) 
    500        
     500 
    501501      ! Apply the masks 
    502502      p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:) 
    503        
     503 
    504504      ! Set missing increments to 0.0 rather than 1e+20 
    505505      ! to allow for differences in masks 
     
    558558         cchl_p_bkg(:,:,:) = 0.0 
    559559#endif 
    560           
     560 
    561561         !-------------------------------------------------------------------- 
    562562         ! Read background variables for phytoplankton assimilation 
     
    578578         CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 
    579579#endif 
    580           
     580 
    581581         IF ( ln_phytobal ) THEN 
    582582 
     
    628628 
    629629         CALL iom_close( inum ) 
    630           
     630 
    631631         DO jt = 1, jptra 
    632632            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    633633         END DO 
    634        
     634 
    635635      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN 
    636636 
     
    641641 
    642642         CALL iom_open( c_asmbkg, inum ) 
    643           
     643 
    644644#if defined key_hadocc 
    645645         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
     
    652652 
    653653         CALL iom_close( inum ) 
    654           
     654 
    655655         DO jt = 1, jptra 
    656656            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    657657         END DO 
    658658         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 
    659        
     659 
    660660      ENDIF 
    661661#else 
     
    681681      !! 
    682682      !! ** Action  : 
    683       !!                    
     683      !! 
    684684      !! References : 
    685685      !! 
     
    695695      REAL(wp)          :: zdate         ! Date 
    696696      !!------------------------------------------------------------------------ 
    697       
     697 
    698698      ! Set things up 
    699699      zdate = REAL( ndastp ) 
     
    706706            &                    TRIM( c_asmbal ) // ' at timestep = ', kt 
    707707 
    708          ! Define the output file        
     708         ! Define the output file 
    709709         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib) 
    710710 
     
    812812            &                    TRIM( c_asmbal ) // ' at timestep = ', kt 
    813813      ENDIF 
    814                                   
     814 
    815815   END SUBROUTINE asm_bgc_bal_wri 
    816816 
     
    893893      !! ** Action  :   return non-log increments 
    894894      !! 
    895       !! References :    
     895      !! References : 
    896896      !!------------------------------------------------------------------------ 
    897897      !! 
     
    970970      !!------------------------------------------------------------------------ 
    971971      !!                    ***  ROUTINE phyto2d_asm_inc  *** 
    972       !!           
     972      !! 
    973973      !! ** Purpose : Apply the chlorophyll assimilation increments. 
    974974      !! 
     
    977977      !!              Direct initialization or Incremental Analysis Updating. 
    978978      !! 
    979       !! ** Action  :  
     979      !! ** Action  : 
    980980      !!------------------------------------------------------------------------ 
    981981      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    10081008      REAL(wp), DIMENSION(jpi,jpj,1) :: zphyt_avg_bkg  ! Local phyt_avg_bkg 
    10091009      !!------------------------------------------------------------------------ 
    1010        
     1010 
    10111011      IF ( kt <= nit000 ) THEN 
    10121012 
     
    10141014         ! Remember that two sets of non-log increments should not be 
    10151015         ! expected to be in the same ratio as their log equivalents 
    1016           
     1016 
    10171017         ! Total chlorophyll 
    10181018         IF ( ln_slchltotinc ) THEN 
     
    11741174 
    11751175            IF(lwp) THEN 
    1176                WRITE(numout,*)  
     1176               WRITE(numout,*) 
    11771177               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', & 
    11781178                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    12051205         ENDIF 
    12061206 
    1207       ELSEIF ( ll_asmdin ) THEN  
     1207      ELSEIF ( ll_asmdin ) THEN 
    12081208 
    12091209         !-------------------------------------------------------------------- 
    12101210         ! Direct Initialization 
    12111211         !-------------------------------------------------------------------- 
    1212           
     1212 
    12131213         IF ( kt == nitdin_r ) THEN 
    12141214 
     
    12341234            END WHERE 
    12351235#endif 
    1236   
     1236 
    12371237            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    12381238            ! which is called at end of model run 
     
    12501250      !!------------------------------------------------------------------------ 
    12511251      !!                    ***  ROUTINE phyto3d_asm_inc  *** 
    1252       !!           
     1252      !! 
    12531253      !! ** Purpose : Apply the profile chlorophyll assimilation increments. 
    12541254      !! 
     
    12561256      !!              Direct initialization or Incremental Analysis Updating. 
    12571257      !! 
    1258       !! ** Action  :  
     1258      !! ** Action  : 
    12591259      !!------------------------------------------------------------------------ 
    12601260      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    13401340 
    13411341            IF(lwp) THEN 
    1342                WRITE(numout,*)  
     1342               WRITE(numout,*) 
    13431343               WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', & 
    13441344                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    13711371         ENDIF 
    13721372 
    1373       ELSEIF ( ll_asmdin ) THEN  
     1373      ELSEIF ( ll_asmdin ) THEN 
    13741374 
    13751375         !-------------------------------------------------------------------- 
    13761376         ! Direct Initialization 
    13771377         !-------------------------------------------------------------------- 
    1378           
     1378 
    13791379         IF ( kt == nitdin_r ) THEN 
    13801380 
     
    14001400            END WHERE 
    14011401#endif 
    1402   
     1402 
    14031403            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    14041404            ! which is called at end of model run 
     
    14171417      !!------------------------------------------------------------------------ 
    14181418      !!                    ***  ROUTINE pco2_asm_inc  *** 
    1419       !!           
     1419      !! 
    14201420      !! ** Purpose : Apply the pco2/fco2 assimilation increments. 
    14211421      !! 
     
    14241424      !!              Direct initialization or Incremental Analysis Updating. 
    14251425      !! 
    1426       !! ** Action  :  
     1426      !! ** Action  : 
    14271427      !!------------------------------------------------------------------------ 
    14281428      INTEGER, INTENT(IN)                   :: kt           ! Current time step 
     
    15861586               jkmax = jpk-1 
    15871587               DO jk = jpk-1, 1, -1 
     1588#if defined key_vvl 
    15881589                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
    15891590                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
     
    15911592                     jkmax = jk 
    15921593                  ENDIF 
     1594#else 
     1595                  IF ( ( zmld(ji,jj) >  gdepw_0(ji,jj,jk)   ) .AND. & 
     1596                     & ( zmld(ji,jj) <= gdepw_0(ji,jj,jk+1) ) ) THEN 
     1597                     zmld(ji,jj) = gdepw_0(ji,jj,jk+1) 
     1598                     jkmax = jk 
     1599                  ENDIF 
     1600#endif 
    15931601               END DO 
    15941602               ! 
     
    16171625 
    16181626            IF(lwp) THEN 
    1619                WRITE(numout,*)  
     1627               WRITE(numout,*) 
    16201628               IF ( ln_spco2inc ) THEN 
    16211629                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 
     
    16541662         ENDIF 
    16551663 
    1656       ELSEIF ( ll_asmdin ) THEN  
     1664      ELSEIF ( ll_asmdin ) THEN 
    16571665 
    16581666         !-------------------------------------------------------------------- 
    16591667         ! Direct Initialization 
    16601668         !-------------------------------------------------------------------- 
    1661           
     1669 
    16621670         IF ( kt == nitdin_r ) THEN 
    16631671 
     
    16831691            END WHERE 
    16841692#endif 
    1685   
     1693 
    16861694            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    16871695            ! which is called at end of model run 
     
    17001708      !!------------------------------------------------------------------------ 
    17011709      !!                    ***  ROUTINE ph_asm_inc  *** 
    1702       !!           
     1710      !! 
    17031711      !! ** Purpose : Apply the pH assimilation increments. 
    17041712      !! 
     
    17071715      !!              Direct initialization or Incremental Analysis Updating. 
    17081716      !! 
    1709       !! ** Action  :  
     1717      !! ** Action  : 
    17101718      !!------------------------------------------------------------------------ 
    17111719      INTEGER,                          INTENT(IN) :: kt        ! Current time step 
     
    17171725      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments 
    17181726      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments 
    1719        
     1727 
    17201728      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation 
    17211729      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation 
     
    17781786            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) 
    17791787         ENDIF 
    1780           
     1788 
    17811789         ! Account for pCO2 balancing if required 
    17821790         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN 
     
    17841792            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) 
    17851793         ENDIF 
    1786           
     1794 
    17871795         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk 
    17881796         ! This requires three calls to the MOCSY carbonate package 
     
    18491857                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic 
    18501858                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk 
    1851                       
     1859 
    18521860                  ENDIF 
    1853                    
     1861 
    18541862               END DO 
    18551863            END DO 
     
    18571865 
    18581866      ENDIF 
    1859        
     1867 
    18601868      IF ( ll_asmiau ) THEN 
    18611869 
     
    18711879 
    18721880            IF(lwp) THEN 
    1873                WRITE(numout,*)  
     1881               WRITE(numout,*) 
    18741882               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & 
    18751883                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    18931901         ENDIF 
    18941902 
    1895       ELSEIF ( ll_asmdin ) THEN  
     1903      ELSEIF ( ll_asmdin ) THEN 
    18961904 
    18971905         !-------------------------------------------------------------------- 
    18981906         ! Direct Initialization 
    18991907         !-------------------------------------------------------------------- 
    1900           
     1908 
    19011909         IF ( kt == nitdin_r ) THEN 
    19021910 
     
    19131921               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
    19141922            END WHERE 
    1915   
     1923 
    19161924            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    19171925            ! which is called at end of model run 
     
    19191927         ! 
    19201928      ENDIF 
    1921 #endif       
     1929#endif 
    19221930      ! 
    19231931   END SUBROUTINE ph_asm_inc 
     
    19301938      !!---------------------------------------------------------------------- 
    19311939      !!                    ***  ROUTINE dyn_asm_inc  *** 
    1932       !!           
     1940      !! 
    19331941      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments. 
    19341942      !! 
    19351943      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    19361944      !! 
    1937       !! ** Action  :  
     1945      !! ** Action  : 
    19381946      !!---------------------------------------------------------------------- 
    19391947      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    20822090 
    20832091            IF(lwp) THEN 
    2084                WRITE(numout,*)  
     2092               WRITE(numout,*) 
    20852093               WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', & 
    20862094                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    21702178#endif 
    21712179            ENDIF 
    2172             
     2180 
    21732181            IF ( kt == nitiaufin_r ) THEN 
    21742182               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) 
     
    21812189         ENDIF 
    21822190 
    2183       ELSEIF ( ll_asmdin ) THEN  
     2191      ELSEIF ( ll_asmdin ) THEN 
    21842192 
    21852193         !-------------------------------------------------------------------- 
    21862194         ! Direct Initialization 
    21872195         !-------------------------------------------------------------------- 
    2188           
     2196 
    21892197         IF ( kt == nitdin_r ) THEN 
    21902198 
     
    22782286#endif 
    22792287            ENDIF 
    2280   
     2288 
    22812289            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) 
    22822290            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc ) 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r8400 r13355  
    1 MODULE sshwzv    
     1MODULE sshwzv 
    22   !!============================================================================== 
    33   !!                       ***  MODULE  sshwzv  *** 
     
    55   !!============================================================================== 
    66   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code 
    7    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA  
     7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA 
    88   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
     
    1717   !!---------------------------------------------------------------------- 
    1818   USE oce             ! ocean dynamics and tracers variables 
    19    USE dom_oce         ! ocean space and time domain variables  
     19   USE dom_oce         ! ocean space and time domain variables 
    2020   USE sbc_oce         ! surface boundary condition: ocean 
    2121   USE domvvl          ! Variable volume 
     
    2828   USE lib_mpp         ! MPP library 
    2929   USE bdy_oce 
    30    USE bdy_par          
     30   USE bdy_par 
    3131   USE bdydyn2d        ! bdy_ssh routine 
    3232#if defined key_agrif 
    3333   USE agrif_opa_interp 
    3434#endif 
    35 #if defined key_asminc    
     35#if defined key_asminc 
    3636   USE asminc          ! Assimilation increment 
    3737#endif 
     
    5959      !!---------------------------------------------------------------------- 
    6060      !!                ***  ROUTINE ssh_nxt  *** 
    61       !!                    
     61      !! 
    6262      !! ** Purpose :   compute the after ssh (ssha) 
    6363      !! 
     
    7373      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv 
    7474      INTEGER, INTENT(in) ::   kt                      ! time step 
    75       !  
     75      ! 
    7676      INTEGER             ::   jk                      ! dummy loop indices 
    7777      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
     
    8080      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
    8181      ! 
    82       CALL wrk_alloc( jpi, jpj, zhdiv )  
     82      CALL wrk_alloc( jpi, jpj, zhdiv ) 
    8383      ! 
    8484      IF( kt == nit000 ) THEN 
     
    102102#if defined key_vvl 
    103103! Don't directly adjust ssh but change hdivn at all levels instead 
    104 ! In trasbc also add in the heat and salt content associated with these changes at each level   
    105         DO jk = 1, jpkm1                                  
    106                  hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk)  
     104! In trasbc also add in the heat and salt content associated with these changes at each level 
     105        DO jk = 1, jpkm1 
     106                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
    107107        END DO 
    108       ENDIF 
    109 #endif 
     108#endif 
     109      ENDIF 
    110110#endif 
    111111 
     
    121121      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 
    122122      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    123       !  
     123      ! 
    124124      z1_rau0 = 0.5_wp * r1_rau0 
    125125      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     
    146146      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    147147      ! 
    148       CALL wrk_dealloc( jpi, jpj, zhdiv )  
     148      CALL wrk_dealloc( jpi, jpj, zhdiv ) 
    149149      ! 
    150150      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
     
    152152   END SUBROUTINE ssh_nxt 
    153153 
    154     
     154 
    155155   SUBROUTINE wzv( kt ) 
    156156      !!---------------------------------------------------------------------- 
    157157      !!                ***  ROUTINE wzv  *** 
    158       !!                    
     158      !! 
    159159      !! ** Purpose :   compute the now vertical velocity 
    160160      !! 
    161       !! ** Method  : - Using the incompressibility hypothesis, the vertical  
    162       !!      velocity is computed by integrating the horizontal divergence   
     161      !! ** Method  : - Using the incompressibility hypothesis, the vertical 
     162      !!      velocity is computed by integrating the horizontal divergence 
    163163      !!      from the bottom to the surface minus the scale factor evolution. 
    164164      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
     
    176176      REAL(wp)            ::   z1_2dt       ! local scalars 
    177177      !!---------------------------------------------------------------------- 
    178        
     178 
    179179      IF( nn_timing == 1 )  CALL timing_start('wzv') 
    180180      ! 
     
    195195      ! 
    196196      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases 
    197          CALL wrk_alloc( jpi, jpj, jpk, zhdiv )  
     197         CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 
    198198         ! 
    199199         DO jk = 1, jpkm1 
     
    215215         END DO 
    216216         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
    217          CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )  
     217         CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 
    218218      ELSE   ! z_star and linear free surface cases 
    219219         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     
    241241      !!                    ***  ROUTINE ssh_nxt  *** 
    242242      !! 
    243       !! ** Purpose :   achieve the sea surface  height time stepping by  
     243      !! ** Purpose :   achieve the sea surface  height time stepping by 
    244244      !!              applying Asselin time filter and swapping the arrays 
    245       !!              ssha  already computed in ssh_nxt   
     245      !!              ssha  already computed in ssh_nxt 
    246246      !! 
    247247      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r11442 r13355  
    2121   USE lib_mpp        ! MPP library 
    2222   USE wrk_nemo       ! work arrays 
    23    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2424   USE stopack 
    2525 
     
    3131 
    3232   INTEGER  ::   albd_init = 0      !: control flag for initialization 
    33    
     33 
    3434   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude 
    3535   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     
    3737   REAL(wp) ::   c2       = 0.10    !  "        " 
    3838   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
    39   
     39 
    4040   !                             !!* namelist namsbc_alb 
    4141   INTEGER  ::   nn_ice_alb 
     
    5252      !!---------------------------------------------------------------------- 
    5353      !!               ***  ROUTINE albedo_ice  *** 
    54       !!           
    55       !! ** Purpose :   Computation of the albedo of the snow/ice system  
    56       !!        
     54      !! 
     55      !! ** Purpose :   Computation of the albedo of the snow/ice system 
     56      !! 
    5757      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
    5858      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
     
    7373      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
    7474      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
    75       !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger  
     75      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
    7676      !!                     under melting conditions than under freezing conditions 
    77       !!                  3) the evolution of ice albedo as a function of ice thickness shows   
     77      !!                  3) the evolution of ice albedo as a function of ice thickness shows 
    7878      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 
    7979      !! 
    8080      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 
    8181      !!                Brandt et al. 2005, J. Climate, vol 18 
    82       !!                Grenfell & Perovich 2004, JGR, vol 109  
     82      !!                Grenfell & Perovich 2004, JGR, vol 109 
    8383      !!---------------------------------------------------------------------- 
    8484      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin) 
     
    9797 
    9898      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    99        
     99 
    100100      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    101101 
    102       IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    103  
    104        
     102      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
     103 
     104 
    105105      SELECT CASE ( nn_ice_alb ) 
    106106 
     
    109109      !------------------------------------------ 
    110110      CASE( 0 ) 
    111         
     111 
    112112         ralb_sf = 0.80       ! dry snow 
    113113         ralb_sm = 0.65       ! melting snow 
    114114         ralb_if = 0.72       ! bare frozen ice 
    115          ralb_im = rn_albice  ! bare puddled ice  
    116           
     115         ralb_im = rn_albice  ! bare puddled ice 
     116 
    117117         !  Computation of ice albedo (free of snow) 
    118118         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
    119119         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if 
    120120         END  WHERE 
    121        
     121 
    122122         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    123123         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
     
    127127         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
    128128         END WHERE 
    129       
     129 
    130130         DO jl = 1, ijpl 
    131131            DO jj = 1, jpj 
     
    133133                  ! freezing snow 
    134134                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
    135                   !                                        !  freezing snow         
     135                  !                                        !  freezing snow 
    136136                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    137137                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
    138138                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   & 
    139                      &        +         zswitch   * ralb_sf   
     139                     &        +         zswitch   * ralb_sf 
    140140 
    141141                  ! melting snow 
     
    143143                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
    144144                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   & 
    145                       &     +         zswitch   *   ralb_sm  
     145                      &     +         zswitch   *   ralb_sm 
    146146                  ! 
    147147                  ! snow albedo 
    148                   zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     148                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 
    149149                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    150                 
     150 
    151151                  ! Ice/snow albedo 
    152152                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     
    155155               END DO 
    156156            END DO 
    157              
     157 
     158#if defined key_traldf_c2d || key_traldf_c3d 
    158159            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
    159160               & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 
    160                          
     161#else 
     162            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
     163               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 
     164                                'key_traldf_c2d or key_traldf_c3d') 
     165#endif 
    161166         END DO 
    162167 
     
    166171      !  New parameterization (2016) 
    167172      !------------------------------------------ 
    168       CASE( 1 )  
     173      CASE( 1 ) 
    169174 
    170175         ralb_im = rn_albice  ! bare puddled ice 
     
    181186!         ralb_sm = 0.82      ! melting snow 
    182187!         ralb_if = 0.54      ! bare frozen ice 
    183 !  
     188! 
    184189         !  Computation of ice albedo (free of snow) 
    185          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     190         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
    186191         z1_c2 = 1. / 0.05 
    187192         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im 
    188193         ELSE WHERE                                              ;   zalb = ralb_if 
    189194         END  WHERE 
    190           
     195 
    191196         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    192197         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
     
    205210 
    206211                   ! snow albedo 
    207                   zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     212                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 
    208213                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    209214 
    210                   ! Ice/snow albedo    
     215                  ! Ice/snow albedo 
    211216                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    212217                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
     
    214219              END DO 
    215220            END DO 
    216              
     221 
     222#if defined key_traldf_c2d || key_traldf_c3d 
    217223            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
    218224               & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 
    219              
     225#else 
     226            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
     227               & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 
     228                                'key_traldf_c2d or key_traldf_c3d') 
     229#endif 
    220230         END DO 
    221231         ! Effect of the clouds (2d order polynomial) 
    222          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );  
     232         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
    223233 
    224234      END SELECT 
    225        
     235 
    226236      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    227237      ! 
     
    232242      !!---------------------------------------------------------------------- 
    233243      !!               ***  ROUTINE albedo_oce  *** 
    234       !!  
     244      !! 
    235245      !! ** Purpose :   Computation of the albedo of the ocean 
    236246      !!---------------------------------------------------------------------- 
     
    238248      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky 
    239249      !! 
    240       REAL(wp) :: zcoef  
     250      REAL(wp) :: zcoef 
    241251      !!---------------------------------------------------------------------- 
    242252      ! 
    243253      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982 
    244       pa_oce_cs(:,:) = zcoef  
     254      pa_oce_cs(:,:) = zcoef 
    245255      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972 
    246256      ! 
     
    257267      !!---------------------------------------------------------------------- 
    258268      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    259       NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
     269      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
    260270      !!---------------------------------------------------------------------- 
    261271      ! 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r11442 r13355  
    161161         !                                      ! ====================== ! 
    162162         ! 
    163          rn_sfac = 1._wp       ! Default to one if missing from namelist  
     163         rn_sfac = 1._wp       ! Default to one if missing from namelist 
    164164         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
    165165         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) 
     
    205205      ENDIF 
    206206 
     207#if defined key_traldf_c2d || key_traldf_c3d 
    207208      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 
    208209         rn_vfac0(:,:) = rn_vfac 
    209210         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 
    210211      ENDIF 
     212#else 
     213      IF ( ln_stopack .AND. nn_spp_relw > 0 ) & 
     214         & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// & 
     215                          'key_traldf_c2d or key_traldf_c3d') 
     216#endif 
    211217 
    212218      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
     
    217223#if defined key_cice 
    218224      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    219          qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)  
     225         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
    220226         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    221227         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
    222228         ENDIF 
    223          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)          
     229         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    224230         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
    225231         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
     
    231237      ! 
    232238   END SUBROUTINE sbc_blk_core 
    233     
    234     
     239 
     240 
    235241   SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 
    236242      !!--------------------------------------------------------------------- 
     
    242248      !! ** Method  :   CORE bulk formulea for the ocean using atmospheric 
    243249      !!      fields read in sbc_read 
    244       !!  
     250      !! 
    245251      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
    246252      !!              - vtau    : j-component of the stress at V-point  (N/m2) 
     
    280286      ! local scalars ( place there for vector optimisation purposes) 
    281287      zcoef_qsatw = 0.98 * 640380. / rhoa 
    282        
     288 
    283289      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    284290 
     
    288294 
    289295      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    290       zwnd_i(:,:) = 0.e0   
     296      zwnd_i(:,:) = 0.e0 
    291297      zwnd_j(:,:) = 0.e0 
    292298#if defined key_cyclone 
     
    332338      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
    333339         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    334      
     340 
    335341      ! ... tau module, i and j component 
    336342      DO jj = 1, jpj 
     
    363369      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    364370 
    365      
     371 
    366372      !  Turbulent fluxes over ocean 
    367373      ! ----------------------------- 
     
    388394         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ') 
    389395      ENDIF 
    390         
     396 
    391397      ! ----------------------------------------------------------------------------- ! 
    392398      !     III    Total FLUXES                                                       ! 
     
    396402         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    397403      ! 
    398       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
     404      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    399405         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    400406         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
     
    437443      ! 
    438444   END SUBROUTINE blk_oce_core 
    439   
    440     
     445 
     446 
    441447#if defined key_lim2 || defined key_lim3 
    442448   SUBROUTINE blk_ice_core_tau 
     
    531537 
    532538      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
    533        
     539 
    534540   END SUBROUTINE blk_ice_core_tau 
    535541 
     
    544550      !!                between atmosphere and sea-ice using CORE bulk 
    545551      !!                formulea, ice variables and read atmmospheric fields. 
    546       !!  
     552      !! 
    547553      !! caution : the net upward water flux has with mm/day unit 
    548554      !!--------------------------------------------------------------------- 
     
    564570      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
    565571      ! 
    566       CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     572      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    567573 
    568574      ! local scalars ( place there for vector optimisation purposes) 
     
    587593               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    588594               ! lw sensitivity 
    589                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     595               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
    590596 
    591597               ! ----------------------------! 
     
    597603               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    598604               ! Latent Heat 
    599                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     605               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   & 
    600606                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    601607              ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    628634 
    629635#if defined  key_lim3 
    630       CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     636      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
    631637 
    632638      ! --- evaporation --- ! 
     
    638644      ! --- evaporation minus precipitation --- ! 
    639645      zsnw(:,:) = 0._wp 
    640       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     646      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
    641647      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    642648      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     
    661667      DO jl = 1, jpl 
    662668         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
    663                                    ! But we do not have Tice => consider it at 0 degC => evap=0  
     669                                   ! But we do not have Tice => consider it at 0 degC => evap=0 
    664670      END DO 
    665671 
    666       CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     672      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
    667673#endif 
    668674 
     
    688694      ! 
    689695      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
    690        
     696 
    691697   END SUBROUTINE blk_ice_core_flx 
    692698#endif 
     
    701707      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    702708      !! 
    703       !! ** Method : Monin Obukhov Similarity Theory  
     709      !! ** Method : Monin Obukhov Similarity Theory 
    704710      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    705711      !! 
     
    711717      !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    712718      !!       rather than temperature difference only... 
    713       !!    - added function "cd_neutral_10m" that uses the improved parametrization of  
     719      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    714720      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    715721      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
     
    746752 
    747753      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z') 
    748      
     754 
    749755      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
    750756      CALL wrk_alloc( jpi,jpj, zeta_u, stab ) 
     
    758764      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    759765 
    760       !! First guess of stability:  
     766      !! First guess of stability: 
    761767      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 
    762768      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     
    772778      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    773779      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    774      
     780 
    775781      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    776782      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10 
     
    783789         ! 
    784790         ztmp1 = T_zu - sst   ! Updating air/sea differences 
    785          ztmp2 = q_zu - q_sat  
     791         ztmp2 = q_zu - q_sat 
    786792 
    787793         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    788794         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta* 
    789795         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q* 
    790         
     796 
    791797         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 
    792798 
    793799         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    794          ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu)  
     800         ztmp0 =  (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 
    795801         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu) 
    796802 
     
    799805         zpsi_h_u = psi_h( zeta_u ) 
    800806         zpsi_m_u = psi_m( zeta_u ) 
    801         
     807 
    802808         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
    803809         IF ( .NOT. l_zt_equal_zu ) THEN 
     
    808814            q_zu = max(0., q_zu) 
    809815         END IF 
    810         
     816 
    811817         IF( ln_cdgw ) THEN      ! surface wave case 
    812             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
     818            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
    813819            Cd      = sqrt_Cd * sqrt_Cd 
    814820         ELSE 
     
    820826           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10 
    821827           sqrt_Cd_n10 = sqrt(ztmp0) 
    822         
     828 
    823829           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b) 
    824830           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
     
    827833           !! Update of transfer coefficients: 
    828834           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a) 
    829            Cd      = ztmp0 / ( ztmp1*ztmp1 )    
     835           Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    830836           sqrt_Cd = SQRT( Cd ) 
    831837         ENDIF 
     
    833839         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    834840         ztmp2 = sqrt_Cd / sqrt_Cd_n10 
    835          ztmp1 = 1. + Ch_n10*ztmp0                
     841         ztmp1 = 1. + Ch_n10*ztmp0 
    836842         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    837843         ! 
    838          ztmp1 = 1. + Ce_n10*ztmp0                
     844         ztmp1 = 1. + Ce_n10*ztmp0 
    839845         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    840846         ! 
     
    854860   FUNCTION cd_neutral_10m( zw10 ) 
    855861      !!---------------------------------------------------------------------- 
    856       !! Estimate of the neutral drag coefficient at 10m as a function  
     862      !! Estimate of the neutral drag coefficient at 10m as a function 
    857863      !! of neutral wind  speed at 10m 
    858864      !! 
     
    860866      !! 
    861867      !! Author: L. Brodeau, june 2014 
    862       !!----------------------------------------------------------------------     
     868      !!---------------------------------------------------------------------- 
    863869      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s) 
    864870      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m 
    865871      ! 
    866872      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33 
    867       !!----------------------------------------------------------------------     
     873      !!---------------------------------------------------------------------- 
    868874      ! 
    869875      CALL wrk_alloc( jpi,jpj, rgt33 ) 
    870876      ! 
    871877      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    872       rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
     878      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1 
    873879      cd_neutral_10m = 1.e-3 * ( & 
    874880         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r11442 r13355  
    2424   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2525   USE timing         ! Timing 
    26    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     26   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2727   USE stopack 
    2828   USE wrk_nemo       ! Memory Allocation 
     
    4242   REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS 
    4343   REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS 
    44    LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term  
     44   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
    4545   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day] 
    4646 
     
    101101               rn_dqdt_s(:,:) = rn_dqdt 
    102102 
    103                IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 
    104                   & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 
     103#if defined key_traldf_c2d || key_traldf_c3d 
     104            IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 
     105               & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 
     106#else 
     107            IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) & 
     108               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 
     109                                'key_traldf_c2d or key_traldf_c3d') 
     110#endif 
     111 
    105112               DO jj = 1, jpj 
    106113                  DO ji = 1, jpi 
     
    117124               CALL wrk_alloc( jpi, jpj, zsrp) 
    118125               zsrp(:,:) = rn_deds 
     126#if defined key_traldf_c2d || key_traldf_c3d 
    119127               IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 
    120128                  & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 
     129#else 
     130            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
     131               & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 
     132                                'key_traldf_c2d or key_traldf_c3d') 
     133#endif 
     134 
     135 
    121136!CDIR COLLAPSE 
    122137               DO jj = 1, jpj 
    123138                  DO ji = 1, jpi 
    124139                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    125                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )  
     140                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
    126141                     sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux 
    127142                     erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only) 
     
    134149               CALL wrk_alloc( jpi, jpj, zsrp) 
    135150               zsrp(:,:) = rn_deds 
     151#if defined key_traldf_c2d || key_traldf_c3d 
    136152               IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 
    137153                  & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 
    138                zerp_bnd = rn_sssr_bnd / rday                          !       -              -     
     154#else 
     155               IF ( ln_stopack .AND. nn_spp_dedt > 0 ) & 
     156                  & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 
     157                                   'key_traldf_c2d or key_traldf_c3d') 
     158#endif 
     159               zerp_bnd = rn_sssr_bnd / rday                          !       -              - 
    139160!CDIR COLLAPSE 
    140161               DO jj = 1, jpj 
    141                   DO ji = 1, jpi                             
     162                  DO ji = 1, jpi 
    142163                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    143164                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
     
    161182   END SUBROUTINE sbc_ssr 
    162183 
    163   
     184 
    164185   SUBROUTINE sbc_ssr_init 
    165186      !!--------------------------------------------------------------------- 
     
    184205      !!---------------------------------------------------------------------- 
    185206      ! 
    186   
    187       REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :  
     207 
     208      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist : 
    188209      READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901) 
    189210901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp ) 
     
    240261      ENDIF 
    241262      ! 
    242       !                            !* Initialize qrp and erp if no restoring  
     263      !                            !* Initialize qrp and erp if no restoring 
    243264      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp 
    244265      IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp 
    245266      ! 
    246267   END SUBROUTINE sbc_ssr_init 
    247        
     268 
    248269   !!====================================================================== 
    249270END MODULE sbcssr 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90

    r12102 r13355  
    4646#define numnam_ref numnam 
    4747#define numnam_cfg numnam 
    48 #define lwm        lwp    
     48#define lwm        lwp 
    4949#define numond     numout 
    5050 
    51 #define wmask      tmask   
     51#define wmask      tmask 
    5252 
    5353#endif 
     
    6363   !!                                          (SPP, SKEB and sea-ice) 
    6464   !!---------------------------------------------------------------------- 
    65    !!    
     65   !! 
    6666   !!   stopack       : Generate stochastic physics perturbations 
    67    !!    
     67   !! 
    6868   !!                   Method 
    6969   !!                   ====== 
     
    7171   !!       - SPPT (Stochastically perturbed parameterization 
    7272   !!         tendencies )scheme for user-selected trends for 
    73    !!            tracers, momentum and sea-ice  
     73   !!            tracers, momentum and sea-ice 
    7474   !!       - SPP (Schastically perturbed parameters) scheme 
    7575   !!         for some (namelist) parameters 
     
    7777   !!         backscatter energy dissipated numerically or 
    7878   !!            through deep convection. 
    79    !!    
    80    !!    
     79   !! 
     80   !! 
    8181   !!                   Acknowledgements: C3S funded ERGO project 
    82    !!    
     82   !! 
    8383   !!---------------------------------------------------------------------- 
    8484   USE par_kind 
    8585   USE timing          ! Timing 
    8686   USE oce             ! ocean dynamics and tracers variables 
    87    USE dom_oce         ! ocean space and time domain variables  
     87   USE dom_oce         ! ocean space and time domain variables 
    8888   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    8989   USE in_out_manager  ! I/O manager 
     
    104104   USE wrk_nemo 
    105105   USE diaptr 
    106    USE zdf_oce          
     106   USE zdf_oce 
    107107   USE phycst 
    108108 
     
    113113   PUBLIC tra_sppt_collect 
    114114   PUBLIC dyn_sppt_collect 
    115    PUBLIC tra_sppt_apply  
    116    PUBLIC dyn_sppt_apply  
     115   PUBLIC tra_sppt_apply 
     116   PUBLIC dyn_sppt_apply 
    117117   PUBLIC stopack_rst 
    118118   PUBLIC stopack_init 
     
    140140   REAL(wp), SAVE :: rn_spp_tau, rn_spp_stdev 
    141141   INTEGER :: skeb_filter_pass, spp_filter_pass 
    142   
     142 
    143143   ! SPPT Logical switches for individual tendencies 
    144144   LOGICAL :: ln_sppt_taumap, ln_stopack_restart, ln_distcoast, & 
    145145   ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, & 
    146146   ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, & 
    147    ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf     
     147   ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf 
    148148   LOGICAL :: & 
    149149   ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad,& 
     
    181181   INTEGER, PARAMETER, PUBLIC :: jk_spp_qsi0   = 8 
    182182   INTEGER, PARAMETER, PUBLIC :: jk_spp_bfr    = 9 
    183    INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd   = 10  
     183   INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd   = 10 
    184184   INTEGER, PARAMETER, PUBLIC :: jk_spp_avt    = 11 
    185185   INTEGER, PARAMETER, PUBLIC :: jk_spp_avm    = 12 
     
    219219   INTEGER, SAVE :: numrep        = 602 
    220220   INTEGER, SAVE :: lkt 
    221     
     221 
    222222   ! Randome generator seed 
    223223   INTEGER, SAVE   :: nn_stopack_seed(4) 
     
    275275   !!---------------------------------------------------------------------- 
    276276   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    277    !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $  
     277   !! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $ 
    278278   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    279279   !!---------------------------------------------------------------------- 
     
    289289      !! 
    290290      !! ** Purpose :   Collect tracer tendencies (additive) 
    291       !!                This function is called by the tendency diagnostics  
     291      !!                This function is called by the tendency diagnostics 
    292292      !!                module 
    293293      !!---------------------------------------------------------------------- 
    294294      INTEGER                   , INTENT(in   ) ::   kt      ! time step 
    295295#endif 
    296       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature  
     296      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
    297297      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity 
    298298      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
     
    354354      !! 
    355355      !! ** Purpose :   Collect momentum tendencies (additive) 
    356       !!                This function is called by the tendency diagnostics  
     356      !!                This function is called by the tendency diagnostics 
    357357      !!                module 
    358358      !!---------------------------------------------------------------------- 
    359359      INTEGER                   , INTENT(in   ) ::   kt      ! time step 
    360360#endif 
    361       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature  
     361      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdx   ! Temperature 
    362362      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptrdy   ! Salinity 
    363363      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index 
     
    469469      !! 
    470470      !! ** Purpose :   Apply collinear perturbation to ice fields 
    471       !!                For specific processes coded in LIM2/LIM3  
     471      !!                For specific processes coded in LIM2/LIM3 
    472472      !!---------------------------------------------------------------------- 
    473473      ! 
     
    529529         zicewrk(:,:,jm) = z5     ; jm=jm+1 
    530530         zicewrk(:,:,jm) = z6     ; jm=jm+1 
    531          zicewrk(:,:,jm) = z7    
     531         zicewrk(:,:,jm) = z7 
    532532       ENDIF 
    533533       IF( kopt .EQ. 3 ) THEN 
     
    601601            CALL lbc_lnk( u_ice, 'U', -1. ) 
    602602            CALL lbc_lnk( v_ice, 'V', -1. ) 
    603 #endif          
     603#endif 
    604604#if defined key_lim2   &&   defined key_lim2_vp 
    605605            CALL lbc_lnk( u_ice(:,1:jpj), 'I', -1. ) 
    606606            CALL lbc_lnk( v_ice(:,1:jpj), 'I', -1. ) 
    607 #endif          
     607#endif 
    608608          ENDIF 
    609609          DEALLOCATE ( zicewrk ) 
     
    616616      !!                  ***  ROUTINE spp_gen *** 
    617617      !! 
    618       !! ** Purpose :   Perturbing parameters (generic function)   
     618      !! ** Purpose :   Perturbing parameters (generic function) 
    619619      !!                Given a value of standard deviation, the 2D parameter 
    620620      !!                coeff is perturbed following an additive Normal, 
     
    633633   CHARACTER (LEN=99) :: cstrng 
    634634   INTEGER :: jklev 
     635 
     636#if defined key_traldf_c3d || key_traldf_c2d 
    635637 
    636638   CALL wrk_alloc(jpi,jpj,gauss) 
     
    687689       jklev = klev 
    688690     ELSE 
    689        jklev = 0  
     691       jklev = 0 
    690692     ENDIF 
    691693     CALL spp_stats(kt,kspp,jklev,coeff) 
     
    693695 
    694696   CALL wrk_dealloc(jpi,jpj,gauss) 
     697 
     698#else 
     699   CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 
     700#endif 
    695701 
    696702   END SUBROUTINE 
     
    707713   REAL(wp) :: mi,ma 
    708714   CHARACTER(LEN=16) :: cstr = '                ' 
    709    SELECT CASE ( kp )  
    710      CASE( jk_spp_alb   )  
     715   SELECT CASE ( kp ) 
     716     CASE( jk_spp_alb   ) 
    711717       cstr = 'ALBEDO      ' 
    712      CASE( jk_spp_rhg   )  
    713        cstr = 'ICE RHEOLOGY'  
    714      CASE( jk_spp_relw  )  
    715        cstr = 'RELATIVE WND'  
    716      CASE( jk_spp_dqdt  )  
    717        cstr = 'SST RELAXAT.'  
    718      CASE( jk_spp_deds  )  
    719        cstr = 'SSS RELAXAT.'  
    720      CASE( jk_spp_arnf  )  
    721        cstr = 'RIVER MIXING'  
    722      CASE( jk_spp_geot  )  
    723        cstr = 'GEOTHERM.FLX'  
    724      CASE( jk_spp_qsi0  )  
    725        cstr = 'SOLAR EXTIN.'  
    726      CASE( jk_spp_bfr   )  
    727        cstr = 'BOTTOM FRICT'  
    728      CASE( jk_spp_aevd  )  
    729        cstr = 'EDDY VISCDIF'  
    730      CASE( jk_spp_avt   )  
     718     CASE( jk_spp_rhg   ) 
     719       cstr = 'ICE RHEOLOGY' 
     720     CASE( jk_spp_relw  ) 
     721       cstr = 'RELATIVE WND' 
     722     CASE( jk_spp_dqdt  ) 
     723       cstr = 'SST RELAXAT.' 
     724     CASE( jk_spp_deds  ) 
     725       cstr = 'SSS RELAXAT.' 
     726     CASE( jk_spp_arnf  ) 
     727       cstr = 'RIVER MIXING' 
     728     CASE( jk_spp_geot  ) 
     729       cstr = 'GEOTHERM.FLX' 
     730     CASE( jk_spp_qsi0  ) 
     731       cstr = 'SOLAR EXTIN.' 
     732     CASE( jk_spp_bfr   ) 
     733       cstr = 'BOTTOM FRICT' 
     734     CASE( jk_spp_aevd  ) 
     735       cstr = 'EDDY VISCDIF' 
     736     CASE( jk_spp_avt   ) 
    731737       cstr = 'VERT. DIFFUS' 
    732      CASE( jk_spp_avm   )  
     738     CASE( jk_spp_avm   ) 
    733739       cstr = 'VERT. VISCOS' 
    734      CASE( jk_spp_tkelc )  
     740     CASE( jk_spp_tkelc ) 
    735741       cstr = 'TKE LANGMUIR' 
    736      CASE( jk_spp_tkedf )  
    737        cstr = 'TKE RN_EDIFF'  
    738      CASE( jk_spp_tkeds )  
    739        cstr = 'TKE RN_EDISS'  
    740      CASE( jk_spp_tkebb )  
     742     CASE( jk_spp_tkedf ) 
     743       cstr = 'TKE RN_EDIFF' 
     744     CASE( jk_spp_tkeds ) 
     745       cstr = 'TKE RN_EDISS' 
     746     CASE( jk_spp_tkebb ) 
    741747       cstr = 'TKE RN_EBB  ' 
    742      CASE( jk_spp_tkefr )  
     748     CASE( jk_spp_tkefr ) 
    743749       cstr = 'TKE RN_EFR  ' 
    744      CASE( jk_spp_ahtu  )  
     750     CASE( jk_spp_ahtu  ) 
    745751       cstr = 'TRALDF AHTU ' 
    746      CASE( jk_spp_ahtv  )  
     752     CASE( jk_spp_ahtv  ) 
    747753       cstr = 'TRALDF AHTV ' 
    748      CASE( jk_spp_ahtw  )  
     754     CASE( jk_spp_ahtw  ) 
    749755       cstr = 'TRALDF AHTW ' 
    750      CASE( jk_spp_ahtt  )  
     756     CASE( jk_spp_ahtt  ) 
    751757       cstr = 'TRALDF AHTT ' 
    752758     CASE( jk_spp_ahubbl ) 
     
    795801 
    796802   DO jp=1,jk_spp 
    797      SELECT CASE ( jp )  
    798        CASE( jk_spp_alb   )  
     803     SELECT CASE ( jp ) 
     804       CASE( jk_spp_alb   ) 
    799805         cstr = 'ALBEDO      ' 
    800        CASE( jk_spp_rhg   )  
    801          cstr = 'ICE RHEOLOGY'  
    802        CASE( jk_spp_relw  )  
    803          cstr = 'RELATIVE WND'  
    804        CASE( jk_spp_dqdt  )  
    805          cstr = 'SST RELAXAT.'  
    806        CASE( jk_spp_deds  )  
    807          cstr = 'SSS RELAXAT.'  
    808        CASE( jk_spp_arnf  )  
    809          cstr = 'RIVER MIXING'  
    810        CASE( jk_spp_geot  )  
    811          cstr = 'GEOTHERM.FLX'  
    812        CASE( jk_spp_qsi0  )  
    813          cstr = 'SOLAR EXTIN.'  
    814        CASE( jk_spp_bfr   )  
    815          cstr = 'BOTTOM FRICT'  
    816        CASE( jk_spp_aevd  )  
    817          cstr = 'EDDY VISCDIF'  
    818        CASE( jk_spp_avt   )  
     806       CASE( jk_spp_rhg   ) 
     807         cstr = 'ICE RHEOLOGY' 
     808       CASE( jk_spp_relw  ) 
     809         cstr = 'RELATIVE WND' 
     810       CASE( jk_spp_dqdt  ) 
     811         cstr = 'SST RELAXAT.' 
     812       CASE( jk_spp_deds  ) 
     813         cstr = 'SSS RELAXAT.' 
     814       CASE( jk_spp_arnf  ) 
     815         cstr = 'RIVER MIXING' 
     816       CASE( jk_spp_geot  ) 
     817         cstr = 'GEOTHERM.FLX' 
     818       CASE( jk_spp_qsi0  ) 
     819         cstr = 'SOLAR EXTIN.' 
     820       CASE( jk_spp_bfr   ) 
     821         cstr = 'BOTTOM FRICT' 
     822       CASE( jk_spp_aevd  ) 
     823         cstr = 'EDDY VISCDIF' 
     824       CASE( jk_spp_avt   ) 
    819825         cstr = 'VERT. DIFFUS' 
    820        CASE( jk_spp_avm   )  
     826       CASE( jk_spp_avm   ) 
    821827         cstr = 'VERT. VISCOS' 
    822        CASE( jk_spp_tkelc )  
     828       CASE( jk_spp_tkelc ) 
    823829         cstr = 'TKE LANGMUIR' 
    824        CASE( jk_spp_tkedf )  
    825          cstr = 'TKE RN_EDIFF'  
    826        CASE( jk_spp_tkeds )  
    827          cstr = 'TKE RN_EDISS'  
    828        CASE( jk_spp_tkebb )  
     830       CASE( jk_spp_tkedf ) 
     831         cstr = 'TKE RN_EDIFF' 
     832       CASE( jk_spp_tkeds ) 
     833         cstr = 'TKE RN_EDISS' 
     834       CASE( jk_spp_tkebb ) 
    829835         cstr = 'TKE RN_EBB  ' 
    830        CASE( jk_spp_tkefr )  
     836       CASE( jk_spp_tkefr ) 
    831837         cstr = 'TKE RN_EFR  ' 
    832        CASE( jk_spp_ahtu  )  
     838       CASE( jk_spp_ahtu  ) 
    833839         cstr = 'TRALDF AHTU ' 
    834        CASE( jk_spp_ahtv  )  
     840       CASE( jk_spp_ahtv  ) 
    835841         cstr = 'TRALDF AHTV ' 
    836        CASE( jk_spp_ahtw  )  
     842       CASE( jk_spp_ahtw  ) 
    837843         cstr = 'TRALDF AHTW ' 
    838        CASE( jk_spp_ahtt  )  
     844       CASE( jk_spp_ahtt  ) 
    839845         cstr = 'TRALDF AHTT ' 
    840846       CASE( jk_spp_ahubbl ) 
     
    922928   INTEGER :: jk 
    923929 
     930#if defined key_traldf_c3d || key_traldf_c2d 
     931 
    924932   CALL wrk_alloc(jpi,jpj,gauss) 
    925933 
     
    969977   ENDIF 
    970978 
     979#else 
     980   CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 
     981 
     982#endif 
     983 
    971984   CALL wrk_dealloc(jpi,jpj,gauss) 
    972985 
     
    9971010   REAL(wp) :: zsd,xme 
    9981011   INTEGER :: jk 
     1012 
     1013#if defined key_dynldf_c3d || key_dynldf_c2d 
    9991014 
    10001015   CALL wrk_alloc(jpi,jpj,gauss) 
     
    10461061 
    10471062   CALL wrk_dealloc(jpi,jpj,gauss) 
     1063 
     1064#else 
     1065   CALL ctl_stop('key_traldf_c1d is not a valid key for STO') 
     1066#endif 
    10481067 
    10491068   END SUBROUTINE 
     
    11951214      ! Note: sshn should be staggered before being used. 
    11961215      SELECT CASE ( cd_type ) 
    1197                CASE ( 'T' )   
     1216               CASE ( 'T' ) 
    11981217                jk=1 
    11991218                zv = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:)*zts(:,:,1) ) 
     
    12851304      ! Random noise on 2d field 
    12861305      IF ( istep == 1 ) THEN 
    1287          CALL sppt_rand2d( g2d )  
     1306         CALL sppt_rand2d( g2d ) 
    12881307         CALL lbc_lnk( g2d , 'T', 1._wp) 
    12891308         g2d_save = g2d 
     
    12971316         g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev 
    12981317      ENDIF 
    1299     
     1318 
    13001319      ! Laplacian filter and re-normalization 
    13011320      DO jf = 1, nk 
     
    13141333      ENDIF 
    13151334#endif 
    1316     
     1335 
    13171336      ! AR(1) process and array swap 
    13181337      g2d = a*gb + b*g2d 
     
    13601379        ENDDO 
    13611380      ENDIF 
    1362        
     1381 
    13631382      ! Bound 
    13641383      IF( nn_sppt_step_bound .EQ. 2 ) THEN 
     
    14821501 
    14831502#ifdef NEMO_V34 
    1484       REWIND( numnam )             
     1503      REWIND( numnam ) 
    14851504      READ  ( numnam, namstopack ) 
    14861505#else 
    1487       REWIND( numnam_ref )  
     1506      REWIND( numnam_ref ) 
    14881507      READ  ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901) 
    14891508901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp ) 
    14901509 
    1491       REWIND( numnam_cfg )   
     1510      REWIND( numnam_cfg ) 
    14921511      READ  ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 ) 
    14931512902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp ) 
     
    15681587         WRITE(numout,*) 
    15691588         WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     spp_filter_pass:', spp_filter_pass 
    1570          WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev       
     1589         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev 
    15711590         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_spp_tau   :', rn_spp_tau 
    15721591         WRITE(numout,*) 
    1573          WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr   
    1574          WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd    
    1575          WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt  
    1576          WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd    
    1577          WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt  
    1578          WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd    
    1579          WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt   
    1580          WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd    
    1581          WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm   
    1582          WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd    
     1592         WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr 
     1593         WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd 
     1594         WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt 
     1595         WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd 
     1596         WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt 
     1597         WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd 
     1598         WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt 
     1599         WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd 
     1600         WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm 
     1601         WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd 
    15831602         WRITE(numout,*) '       SPP for solar penetration scheme  (only RGB)        nn_spp_qsi0  :', nn_spp_qsi0 
    1584          WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd    
     1603         WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd 
    15851604         WRITE(numout,*) '       SPP for horiz. diffusivity  U                       nn_spp_ahtu  :', nn_spp_ahtu 
    1586          WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd    
     1605         WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd 
    15871606         WRITE(numout,*) '       SPP for horiz. diffusivity  V                       nn_spp_ahtv  :', nn_spp_ahtv 
    1588          WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd    
     1607         WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd 
    15891608         WRITE(numout,*) '       SPP for horiz. diffusivity  W                       nn_spp_ahtw  :', nn_spp_ahtw 
    1590          WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd    
     1609         WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd 
    15911610         WRITE(numout,*) '       SPP for horiz. diffusivity  T                       nn_spp_ahtt  :', nn_spp_ahtt 
    1592          WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd    
     1611         WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd 
    15931612         WRITE(numout,*) '       SPP for horiz. viscosity (1/3)                      nn_spp_ahm1  :', nn_spp_ahm1 
    1594          WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd    
     1613         WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd 
    15951614         WRITE(numout,*) '       SPP for horiz. viscosity (2/4)                      nn_spp_ahm2  :', nn_spp_ahm2 
    1596          WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd    
     1615         WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd 
    15971616         WRITE(numout,*) '       SPP for relative wind factor                        nn_spp_relw  :', nn_spp_relw 
    15981617         WRITE(numout,*) '       (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]' 
    1599          WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd    
     1618         WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd 
    16001619         WRITE(numout,*) '       SPP for mixing close to river mouth                 nn_spp_arnf  :', nn_spp_arnf 
    1601          WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd    
     1620         WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd 
    16021621         WRITE(numout,*) '       SPP for geothermal heating                          nn_spp_geot  :', nn_spp_geot 
    1603          WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd    
     1622         WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd 
    16041623         WRITE(numout,*) '       SPP for enhanced vertical diffusion                 nn_spp_aevd  :', nn_spp_aevd 
    1605          WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd    
     1624         WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd 
    16061625         WRITE(numout,*) '       SPP for TKE rn_lc    Langmuir cell coefficient      nn_spp_tkelc :', nn_spp_tkelc 
    1607          WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd    
     1626         WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd 
    16081627         WRITE(numout,*) '       SPP for TKE rn_ediff Eddy diff. coefficient         nn_spp_tkedf :', nn_spp_tkedf 
    1609          WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd    
     1628         WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd 
    16101629         WRITE(numout,*) '       SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds 
    1611          WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd    
     1630         WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd 
    16121631         WRITE(numout,*) '       SPP for TKE rn_ebb   Surface input of tke           nn_spp_tkebb :', nn_spp_tkebb 
    1613          WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd    
     1632         WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd 
    16141633         WRITE(numout,*) '       SPP for TKE rn_efr   Fraction of srf TKE below ML   nn_spp_tkefr :', nn_spp_tkefr 
    1615          WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd    
     1634         WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd 
    16161635         WRITE(numout,*) '       SPP for BBL U  diffusivity                          nn_spp_ahubbl:', nn_spp_ahubbl 
    16171636         WRITE(numout,*) '                                            STDEV          rn_ahubbl_sd :', rn_ahubbl_sd 
     
    16261645         WRITE(numout,*) 
    16271646         WRITE(numout,*) ' SKEB Perturbation scheme ' 
    1628          WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb     
    1629          WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb     
     1647         WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb 
     1648         WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb 
    16301649         WRITE(numout,*) '       Frequency update for dissipation mask               nn_dcom_freq :', nn_dcom_freq 
    16311650         WRITE(numout,*) '       Numerical dissipation factor (resolut. dependent)   rn_kh        :', rn_kh 
    16321651         WRITE(numout,*) '       Number of passes for spatial filter (AR1 field)     skeb_filter_pass:', skeb_filter_pass 
    1633          WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev       
     1652         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev 
    16341653         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_skeb_tau  :', rn_skeb_tau 
    16351654         WRITE(numout,*) '       Option of convection energy dissipation             nn_dconv     :', nn_dconv 
     
    17521771 
    17531772      ! Find filter attenuation factor 
    1754     
     1773 
    17551774      flt_fac = sppt_fltfac( sppt_filter_pass ) 
    17561775      rdt_sppt = nn_rndm_freq * rn_rdt 
    1757     
     1776 
    17581777      IF( ln_sppt_taumap ) THEN 
    17591778         CALL iom_open ( 'sppt_tau_map', inum ) 
     
    17981817      gauss_b = 0._wp 
    17991818      ! Weigths 
    1800       gauss_w(:)    = 1.0_wp  
     1819      gauss_w(:)    = 1.0_wp 
    18011820      IF( nn_vwei .eq. 1 ) THEN 
    18021821        gauss_w(1)    = 0.0_wp 
     
    18611880      IF(lwp .and. ln_stopack_diags) & 
    18621881      CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    1863     
     1882 
    18641883   END SUBROUTINE stopack_init 
    18651884   ! 
     
    18741893      INTEGER :: id1, jseed 
    18751894      CHARACTER(LEN=10)   ::   clseed='spsd0_0000' 
    1876       INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type  
    1877       INTEGER(KIND=8)     ::   ivals(8)  
     1895      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type 
     1896      INTEGER(KIND=8)     ::   ivals(8) 
    18781897      REAL(wp)            ::   zrseed4(4)           ! RNG seeds in integer type 
    18791898      REAL(wp)            ::   zrseed2d(jpi,jpj) 
     
    19832002      !!--------------------------------------------------------------------- 
    19842003      ! 
    1985       ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,&  
    1986       gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , &  
     2004      ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 
     2005      gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 
    19872006      spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,& 
    19882007      gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),& 
     
    22082227      IF ( ln_dpsiv ) THEN 
    22092228       DO jp=1,jpni-1 
    2210          IF( jpri == jp ) THEN ! SEND TO EAST  
     2229         IF( jpri == jp ) THEN ! SEND TO EAST 
    22112230          zwrk(1:jpj) = dpsiv(jpi-1,:) 
    22122231          tag=2000+narea 
     
    22682287   REAL :: ds,dt,dtot,kh2 
    22692288   INTEGER :: ji,jj,jk 
    2270     
     2289 
    22712290   IF ( mt .eq. nit000 ) THEN 
    22722291        ALLOCATE ( dnum(jpi,jpj,jpk) ) 
     
    22872306          dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e2v(ji,jj) + & 
    22882307               (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e1u(ji,jj) 
    2289    
     2308 
    22902309          dtot = sqrt( ds*ds + dt*dt ) * tmask(ji,jj,jk) 
    22912310          dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(ji,jj)*e2t(ji,jj) 
     
    22932312      ENDDO 
    22942313     ENDDO 
    2295     
     2314 
    22962315     CALL lbc_lnk(dnum,'T',1._wp) 
    22972316 
    22982317   ENDIF 
    22992318 
    2300    END SUBROUTINE  
     2319   END SUBROUTINE 
    23012320 
    23022321   SUBROUTINE skeb_dcon ( mt ) 
     
    23292348 
    23302349           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) & 
    2331               & / ( rau0 * fse3w(ji,jj,jk) )  
     2350              & / ( rau0 * fse3w(ji,jj,jk) ) 
    23322351 
    23332352           dcon(ji,jj,jk) = kc2*zz*e1t(ji,jj)*e2t(ji,jj)*rau0 / fse3w(ji,jj,jk) 
     
    23782397     IF(ln_skeb_own_gauss) THEN 
    23792398       DO jk=1,jpkm1 
    2380          psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2399         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 
    23812400       ENDDO 
    23822401     ELSE 
     
    24072426     IF(ln_skeb_own_gauss) THEN 
    24082427       DO jk=1,jpkm1 
    2409          psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2428         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
    24102429       ENDDO 
    24112430     ELSE 
     
    24402459   IF(ln_skeb_own_gauss) THEN 
    24412460     DO jk=1,jpkm1 
    2442        psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2461       psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
    24432462     ENDDO 
    24442463   ELSE 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r11442 r13355  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   tra_bbc      : update the tracer trend at ocean bottom  
     14   !!   tra_bbc      : update the tracer trend at ocean bottom 
    1515   !!   tra_bbc_init : initialization of geothermal heat flux trend 
    1616   !!---------------------------------------------------------------------- 
     
    1919   USE phycst          ! physical constants 
    2020   USE trd_oce         ! trends: ocean variables 
    21    USE trdtra          ! trends manager: tracers  
     21   USE trdtra          ! trends manager: tracers 
    2222   USE in_out_manager  ! I/O manager 
    2323   USE iom             ! I/O manager 
     
    4444   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend 
    4545   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read) 
    46   
     46 
    4747   !! * Substitutions 
    4848#  include "domzgr_substitute.h90" 
     
    5858      !!                  ***  ROUTINE tra_bbc  *** 
    5959      !! 
    60       !! ** Purpose :   Compute the bottom boundary contition on temperature  
    61       !!              associated with geothermal heating and add it to the  
     60      !! ** Purpose :   Compute the bottom boundary contition on temperature 
     61      !!              associated with geothermal heating and add it to the 
    6262      !!              general trend of temperature equations. 
    6363      !! 
    64       !! ** Method  :   The geothermal heat flux set to its constant value of  
     64      !! ** Method  :   The geothermal heat flux set to its constant value of 
    6565      !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999). 
    6666      !!       The temperature trend associated to this heat flux through the 
     
    9292      !                             !  Add the geothermal heat flux trend on temperature 
    9393 
     94#if defined key_traldf_c2d || key_traldf_c3d 
    9495      IF( ln_stopack .AND. nn_spp_geot > 0) THEN 
    9596          qgh_trd1(:,:) = qgh_trd0(:,:) 
    9697          CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot) 
    9798      ENDIF 
     99#else 
     100      IF ( ln_stopack .AND. nn_spp_geot > 0 ) & 
     101         & CALL ctl_stop( 'tra_bbc: parameter perturbation will only work with '// & 
     102                          'key_traldf_c2d or key_traldf_c3d') 
     103#endif 
     104 
     105 
    98106      DO jj = 2, jpjm1 
    99107         DO ji = 2, jpim1 
     
    144152      CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files 
    145153      ! 
    146       NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir  
     154      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 
    147155      !!---------------------------------------------------------------------- 
    148156 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r11442 r13355  
    3232   USE trdtra         ! trends: active tracers 
    3333   ! 
    34    USE iom            ! IOM library                
     34   USE iom            ! IOM library 
    3535   USE in_out_manager ! I/O manager 
    3636   USE lbclnk         ! ocean lateral boundary conditions 
     
    3838   USE wrk_nemo       ! Memory Allocation 
    3939   USE timing         ! Timing 
    40    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    4141   USE stopack 
    4242 
     
    199199      ! 
    200200      ahu_bbl_1(:,:) = ahu_bbl(:,:) 
     201#if defined key_traldf_c2d || key_traldf_c3d 
    201202      IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN 
    202203          CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl ) 
    203204      ENDIF 
     205#else 
     206      IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) & 
     207         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// & 
     208                          'key_traldf_c2d or key_traldf_c3d') 
     209#endif 
     210 
     211 
    204212      ahv_bbl_1(:,:) = ahv_bbl(:,:) 
     213#if defined key_traldf_c2d || key_traldf_c3d 
    205214      IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN 
    206215          CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl ) 
    207216      ENDIF 
     217#else 
     218      IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) & 
     219         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// & 
     220                          'key_traldf_c2d or key_traldf_c3d') 
     221#endif 
     222 
    208223      ! 
    209224      DO jn = 1, kjpt                                     ! tracer loop 
     
    215230            END DO 
    216231         END DO 
    217          !                
     232         ! 
    218233         DO jj = 2, jpjm1                                    ! Compute the trend 
    219234            DO ji = 2, jpim1 
     
    426441                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
    427442                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
    428                   !                                                          ! 2*masked bottom density gradient  
     443                  !                                                          ! 2*masked bottom density gradient 
    429444                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
    430445                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
     
    585600               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    586601            ENDIF 
    587             !      
     602            ! 
    588603            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
    589604               mgrhv(ji,jj) = INT(  SIGN( 1.e0, & 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r11442 r13355  
    99   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    1010   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate 
    11    !!            3.2  !  2009-04  (G. Madec & NEMO team)  
    12    !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     11   !!            3.2  !  2009-04  (G. Madec & NEMO team) 
     12   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
    1313   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 
    1414   !!---------------------------------------------------------------------- 
     
    4343   !                                 !!* Namelist namtra_qsr: penetrative solar radiation 
    4444   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag 
    45    LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag   
     45   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
    4646   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag 
    4747   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag 
     
    5151   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    5252   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands) 
    53   
     53 
    5454   ! Module variables 
    5555   REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0 
     
    8080      !!      Considering the 2 wavebands case: 
    8181      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
    82       !!         The temperature trend associated with the solar radiation penetration  
     82      !!         The temperature trend associated with the solar radiation penetration 
    8383      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 
    8484      !!         At the bottom, boudary condition for the radiation is no flux : 
     
    8686      !!      in the last ocean level. 
    8787      !!         In z-coordinate case, the computation is only done down to the 
    88       !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients  
     88      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients 
    8989      !!      used for the computation are calculated one for once as they 
    9090      !!      depends on k only. 
     
    106106      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         - 
    107107      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
    108       REAL(wp) ::   zlogc, zlogc2, zlogc3  
     108      REAL(wp) ::   zlogc, zlogc2, zlogc3 
    109109      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
    110110      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d 
     
    113113      IF( nn_timing == 1 )  CALL timing_start('tra_qsr') 
    114114      ! 
    115       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    116       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )  
     115      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
     116      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
    117117      ! 
    118118      IF( kt == nit000 ) THEN 
     
    124124 
    125125      IF( l_trdtra ) THEN      ! Save ta and sa trends 
    126          CALL wrk_alloc( jpi, jpj, jpk, ztrdt )  
     126         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
    127127         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    128128      ENDIF 
     
    150150      !                                        Compute now qsr tracer content field 
    151151      !                                        ************************************ 
    152        
     152 
    153153      !                                           ! ============================================== ! 
    154154      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
     
    159159         !                                        Add to the general trend 
    160160         DO jk = 1, jpkm1 
    161             DO jj = 2, jpjm1  
     161            DO jj = 2, jpjm1 
    162162               DO ji = fs_2, fs_jpim1   ! vector opt. 
    163163                  z1_e3t = zfact / fse3t(ji,jj,jk) 
     
    180180         ENDIF 
    181181         !                                        ! ============================================== ! 
    182       ELSE                                        !  Ocean alone :  
     182      ELSE                                        !  Ocean alone : 
    183183         !                                        ! ============================================== ! 
    184184         ! 
    185          !  
     185         ! 
     186#if defined key_traldf_c2d || key_traldf_c3d 
    186187         IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN 
    187              xsi0r = rn_si0 
    188              CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 
    189              xsi0r = 1.e0 / xsi0r 
    190          ENDIF 
     188            xsi0r = rn_si0 
     189            CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 
     190            xsi0r = 1.e0 / xsi0r 
     191         ENDIF 
     192#else 
     193         IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) & 
     194            & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// & 
     195                             'key_traldf_c2d or key_traldf_c3d') 
     196#endif 
     197 
    191198         !                                               ! ------------------------- ! 
    192199         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration ! 
     
    199206                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step 
    200207                  DO jk = 1, nksr + 1 
    201                      zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)  
     208                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
    202209                  ENDDO 
    203210                  ! 
     
    220227                        zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
    221228                        zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
    222                         zCze    = 1.12  * (zchl)**0.803  
     229                        zCze    = 1.12  * (zchl)**0.803 
    223230                        DO jk = 1, nksr + 1 
    224231                           zpsi = fsdept(ji,jj,jk) / zze 
     
    231238               ELSE                              !* Variable ocean volume but constant chrlorophyll 
    232239                  DO jk = 1, nksr + 1 
    233                      zchl3d(:,:,jk) = 0.05  
     240                     zchl3d(:,:,jk) = 0.05 
    234241                  ENDDO 
    235242               ENDIF 
     
    256263!CDIR NOVERRCHK 
    257264                  DO jj = 1, jpj 
    258 !CDIR NOVERRCHK    
     265!CDIR NOVERRCHK 
    259266                     DO ji = 1, jpi 
    260267                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) 
     
    289296                     END DO 
    290297                  END DO 
    291                   !  
     298                  ! 
    292299                  DO jj = 1, jpj 
    293300                     DO ji = 1, jpi 
     
    296303                        zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
    297304                        zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
    298                         fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)  
     305                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
    299306                     END DO 
    300307                  END DO 
     
    321328               zz0   =        rn_abs   * r1_rau0_rcp 
    322329               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp 
    323                DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m  
     330               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
    324331                  DO jj = 1, jpj 
    325332                     DO ji = 1, jpi 
    326333                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    327334                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    328                         qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )  
     335                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
    329336                     END DO 
    330337                  END DO 
     
    344351                  DO jj = 2, jpjm1 
    345352                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    346                         ! (ISF) no light penetration below the ice shelves          
     353                        ! (ISF) no light penetration below the ice shelves 
    347354                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 
    348355                     END DO 
     
    360367         !                                        Add to the general trend 
    361368         DO jk = 1, nksr 
    362             DO jj = 2, jpjm1  
     369            DO jj = 2, jpjm1 
    363370               DO ji = fs_2, fs_jpim1   ! vector opt. 
    364371                  z1_e3t = zfact / fse3t(ji,jj,jk) 
     
    376383            &                    'at it= ', kt,' date= ', ndastp 
    377384         IF(lwp) WRITE(numout,*) '~~~~' 
    378          IF(nn_timing == 2)  CALL timing_start('iom_rstput')  
     385         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
    379386         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      ) 
    380          CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm  
     387         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
    381388         IF(nn_timing == 2)  CALL timing_stop('iom_rstput') 
    382389         ! 
     
    386393         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    387394         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    388          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )  
     395         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
    389396      ENDIF 
    390397      !                       ! print mean trends (used for debugging) 
    391398      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
    392399      ! 
    393       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    394       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )  
     400      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
     401      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
    395402      ! 
    396403      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr') 
     
    408415      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio 
    409416      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The 
    410       !!      default values correspond to clear water (type I in Jerlov'  
     417      !!      default values correspond to clear water (type I in Jerlov' 
    411418      !!      (1968) classification. 
    412419      !!         called by tra_qsr at the first timestep (nit000) 
     
    435442      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init') 
    436443      ! 
    437       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    438       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     444      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
     445      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
    439446      ! 
    440447 
     
    465472 
    466473      IF( ln_traqsr ) THEN     ! control consistency 
    467          !                       
     474         ! 
    468475         IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN 
    469476            CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 
     
    480487            &              ' 2 bands, 3 RGB bands or bio-model light penetration' ) 
    481488         ! 
    482          IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1  
     489         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
    483490         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2 
    484491         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3 
     
    497504      ENDIF 
    498505      !                          ! ===================================== ! 
    499       IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !   
     506      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
    500507         !                       ! ===================================== ! 
    501508         ! 
     
    539546                  zchl = 0.05                                 ! constant chlorophyll 
    540547                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    541                   zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll  
     548                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll 
    542549                  zekg(:,:) = rkrgb(2,irgb) 
    543550                  zekr(:,:) = rkrgb(3,irgb) 
     
    546553                  ze0(:,:,1) = rn_abs 
    547554                  ze1(:,:,1) = zcoef 
    548                   ze2(:,:,1) = zcoef  
     555                  ze2(:,:,1) = zcoef 
    549556                  ze3(:,:,1) = zcoef 
    550557                  zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask 
    551                 
     558 
    552559                  DO jk = 2, nksr+1 
    553560!CDIR NOVERRCHK 
    554561                     DO jj = 1, jpj 
    555 !CDIR NOVERRCHK    
     562!CDIR NOVERRCHK 
    556563                        DO ji = 1, jpi 
    557564                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) 
     
    566573                        END DO 
    567574                     END DO 
    568                   END DO  
     575                  END DO 
    569576                  ! 
    570577                  DO jk = 1, nksr 
     
    598605                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    599606                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    600                         etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)  
     607                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
    601608                     END DO 
    602609                  END DO 
     
    607614         ENDIF 
    608615         !                       ! ===================================== ! 
    609       ELSE                       !        No light penetration           !                    
     616      ELSE                       !        No light penetration           ! 
    610617         !                       ! ===================================== ! 
    611618         IF(lwp) THEN 
     
    625632      ENDIF 
    626633      ! 
    627       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    628       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     634      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
     635      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
    629636      ! 
    630637      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init') 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r11442 r13355  
    107107         WRITE(numout,*) '~~~~~~~~' 
    108108      ENDIF 
    109       !  
     109      ! 
     110#if defined key_traldf_c2d || key_traldf_c3d 
    110111      IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN 
    111112         bfrcoef2d = bfrcoef2d0 
    112113         CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr) 
    113114      ENDIF 
     115#else 
     116      IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) & 
     117         & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// & 
     118                          'key_traldf_c2d or key_traldf_c3d') 
     119#endif 
    114120      ! 
    115121      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only 
     
    140146               END DO 
    141147            END IF 
    142          !    
     148         ! 
    143149         ELSE 
    144150            zbfrt(:,:) = bfrcoef2d(:,:) 
     
    183189               DO ji = 2, jpim1 
    184190                  ! (ISF) ======================================================================== 
    185                   ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     191                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
    186192                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points) 
    187193                  ! 
     
    363369            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    364370         ENDIF 
    365           
     371 
    366372         IF ( ln_isfcav ) THEN 
    367373            IF(ln_tfr2d) THEN 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r11442 r13355  
    2020   USE zdfkpp          ! KPP vertical mixing 
    2121   USE trd_oce         ! trends: ocean variables 
    22    USE trdtra          ! trends manager: tracers  
     22   USE trdtra          ! trends manager: tracers 
    2323   USE in_out_manager  ! I/O manager 
    2424   USE iom             ! for iom_put 
    2525   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2626   USE timing          ! Timing 
    27    USE stopack          
     27   USE stopack 
    2828 
    2929   IMPLICIT NONE 
     
    4545      !!---------------------------------------------------------------------- 
    4646      !!                  ***  ROUTINE zdf_evd  *** 
    47       !!                    
     47      !! 
    4848      !! ** Purpose :   Local increased the vertical eddy viscosity and diffu- 
    4949      !!      sivity coefficients when a static instability is encountered. 
    5050      !! 
    5151      !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients 
    52       !!      are set to avevd (namelist parameter) if the water column is  
     52      !!      are set to avevd (namelist parameter) if the water column is 
    5353      !!      statically unstable (i.e. if rn2 < -1.e-12 ) 
    5454      !! 
     
    7777      zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
    7878 
     79#if defined key_traldf_c2d || key_traldf_c3d 
    7980      IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN 
    8081         rn_avevd0(:,:) = rn_avevd 
    8182         CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd) 
    8283      ENDIF 
     84#else 
     85      IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) & 
     86         & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// & 
     87                          'key_traldf_c2d or key_traldf_c3d') 
     88#endif 
    8389 
    8490      SELECT CASE ( nn_evdm ) 
     
    8894         zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application 
    8995         ! 
    90          DO jk = 1, jpkm1  
     96         DO jk = 1, jpkm1 
    9197            DO jj = 2, jpj             ! no vector opt. 
    9298               DO ji = 2, jpi 
     
    106112               END DO 
    107113            END DO 
    108          END DO  
     114         END DO 
    109115         CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions 
    110116         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. ) 
     
    113119         CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
    114120         ! 
    115       CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12)  
     121      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
    116122         DO jk = 1, jpkm1 
    117 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!  
     123!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
    118124            DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call) 
    119125               DO ji = 1, jpi 
    120126#if defined key_zdfkpp 
    121127                  ! no evd mixing in the boundary layer with KPP 
    122                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &           
     128                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   & 
    123129#else 
    124130                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   & 
     
    129135         END DO 
    130136         ! 
    131       END SELECT  
     137      END SELECT 
    132138 
    133139      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r11442 r13355  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  zdfgls  *** 
    4    !! Ocean physics:  vertical mixing coefficient computed from the gls  
     4   !! Ocean physics:  vertical mixing coefficient computed from the gls 
    55   !!                 turbulent closure parameterization 
    66   !!====================================================================== 
     
    1616   !!   gls_rst       : read/write gls restart in ocean restart file 
    1717   !!---------------------------------------------------------------------- 
    18    USE oce            ! ocean dynamics and active tracers  
     18   USE oce            ! ocean dynamics and active tracers 
    1919   USE dom_oce        ! ocean space and time domain 
    2020   USE domvvl         ! ocean space and time domain : variable volume layer 
     
    3131   USE iom            ! I/O manager library 
    3232   USE timing         ! Timing 
    33    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     33   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3434   USE stopack 
    3535 
     
    6262   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
    6363   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    64    REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
     64   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
    6565 
    6666   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
    67    REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1  
    68    REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn     
     67   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
     68   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn 
    6969   REAL(wp) ::   rghmin        = -0.28_wp 
    7070   REAL(wp) ::   rgh0          =  0.0329_wp 
     
    7373   REAL(wp) ::   ra2           =  0.74_wp 
    7474   REAL(wp) ::   rb1           = 16.60_wp 
    75    REAL(wp) ::   rb2           = 10.10_wp          
    76    REAL(wp) ::   re2           =  1.33_wp          
     75   REAL(wp) ::   rb2           = 10.10_wp 
     76   REAL(wp) ::   re2           =  1.33_wp 
    7777   REAL(wp) ::   rl1           =  0.107_wp 
    7878   REAL(wp) ::   rl2           =  0.0032_wp 
     
    134134      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments 
    135135      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars 
    136       REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -  
     136      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
    137137      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      - 
    138138      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      - 
     
    140140      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep 
    141141      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar 
    142       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves  
     142      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
    143143      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves) 
    144144      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before 
     
    157157      CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 
    158158      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  ) 
    159        
     159 
    160160      ! Preliminary computing 
    161161 
     
    166166         avm (:,:,:) = avm_k (:,:,:) 
    167167         avmu(:,:,:) = avmu_k(:,:,:) 
    168          avmv(:,:,:) = avmv_k(:,:,:)  
     168         avmv(:,:,:) = avmv_k(:,:,:) 
    169169      ENDIF 
    170170 
    171171      ! Compute surface and bottom friction at T-points 
    172 !CDIR NOVERRCHK           
    173       DO jj = 2, jpjm1           
    174 !CDIR NOVERRCHK          
    175          DO ji = fs_2, fs_jpim1   ! vector opt.          
     172!CDIR NOVERRCHK 
     173      DO jj = 2, jpjm1 
     174!CDIR NOVERRCHK 
     175         DO ji = fs_2, fs_jpim1   ! vector opt. 
    176176            ! 
    177177            ! surface friction 
    178178            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
    179             !    
    180             ! bottom friction (explicit before friction)         
    181             ! Note that we chose here not to bound the friction as in dynbfr)    
    182             ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &          
    183                & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )       
    184             zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &          
    185                & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )       
    186             ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)          
    187          END DO          
    188       END DO     
     179            ! 
     180            ! bottom friction (explicit before friction) 
     181            ! Note that we chose here not to bound the friction as in dynbfr) 
     182            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   & 
     183               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  ) 
     184            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   & 
     185               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  ) 
     186            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
     187         END DO 
     188      END DO 
    189189 
    190190      ! Set surface roughness length 
    191191      SELECT CASE ( nn_z0_met ) 
    192192      ! 
    193       CASE ( 0 )             ! Constant roughness           
     193      CASE ( 0 )             ! Constant roughness 
    194194         zhsro(:,:) = rn_hsro 
    195195      CASE ( 1 )             ! Standard Charnock formula 
     
    227227      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    228228         DO jk = 2, jpkm1 
    229             DO jj = 2, jpjm1  
     229            DO jj = 2, jpjm1 
    230230               DO ji = fs_2, fs_jpim1   ! vector opt. 
    231231                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) 
     
    276276               IF( ln_sigpsi ) THEN 
    277277                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1. 
    278                   zwall_psi(ji,jj,jk) = rsc_psi /   &  
     278                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
    279279                     &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  ) 
    280280               ELSE 
     
    295295               ! diagonal 
    296296               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
    297                   &                       + rdt * zdiss * tmask(ji,jj,jk)  
     297                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
    298298               ! 
    299299               ! right hand side in en 
     
    317317      ! First level 
    318318      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    319       en(:,:,1) = MAX(en(:,:,1), rn_emin)  
     319      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
    320320      z_elem_a(:,:,1) = en(:,:,1) 
    321321      z_elem_c(:,:,1) = 0._wp 
    322322      z_elem_b(:,:,1) = 1._wp 
    323       !  
     323      ! 
    324324      ! One level below 
    325325      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    326326          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    327327      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    328       z_elem_a(:,:,2) = 0._wp  
     328      z_elem_a(:,:,2) = 0._wp 
    329329      z_elem_c(:,:,2) = 0._wp 
    330330      z_elem_b(:,:,2) = 1._wp 
     
    335335      ! Dirichlet conditions at k=1 
    336336      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    337       en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
     337      en(:,:,1)       = MAX(en(:,:,1), rn_emin) 
    338338      z_elem_a(:,:,1) = en(:,:,1) 
    339339      z_elem_c(:,:,1) = 0._wp 
     
    358358      SELECT CASE ( nn_bc_bot ) 
    359359      ! 
    360       CASE ( 0 )             ! Dirichlet  
     360      CASE ( 0 )             ! Dirichlet 
    361361         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 
    362362         !                      ! Balance between the production and the dissipation terms 
     
    378378               z_elem_c(ji,jj,ibotm1) = 0._wp 
    379379               z_elem_b(ji,jj,ibotm1) = 1._wp 
    380                en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )  
     380               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    381381            END DO 
    382382         END DO 
    383383         ! 
    384384      CASE ( 1 )             ! Neumman boundary condition 
    385          !                       
     385         ! 
    386386!CDIR NOVERRCHK 
    387387         DO jj = 2, jpjm1 
     
    429429         END DO 
    430430      END DO 
    431       !                                            ! set the minimum value of tke  
     431      !                                            ! set the minimum value of tke 
    432432      en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
    433433 
     
    471471            DO jj = 2, jpjm1 
    472472               DO ji = fs_2, fs_jpim1   ! vector opt. 
    473                   psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn  
     473                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
    474474               END DO 
    475475            END DO 
     
    490490               ! 
    491491               ! psi / k 
    492                zratio = psi(ji,jj,jk) / eb(ji,jj,jk)  
     492               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
    493493               ! 
    494494               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) 
     
    510510               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term 
    511511               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
    512                !                                                         
     512               ! 
    513513               ! building the matrix 
    514514               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
     
    552552      z_elem_c(:,:,2) = 0._wp 
    553553      z_elem_b(:,:,2) = 1._wp 
    554       !  
     554      ! 
    555555      ! 
    556556      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
     
    576576      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    577577 
    578       !    
     578      ! 
    579579      ! 
    580580      END SELECT 
     
    586586      ! 
    587587      ! 
    588       CASE ( 0 )             ! Dirichlet  
     588      CASE ( 0 )             ! Dirichlet 
    589589         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 
    590590         !                      ! Balance between the production and the dissipation terms 
     
    611611         ! 
    612612      CASE ( 1 )             ! Neumman boundary condition 
    613          !                       
     613         ! 
    614614!CDIR NOVERRCHK 
    615615         DO jj = 2, jpjm1 
     
    693693            DO jj = 2, jpjm1 
    694694               DO ji = fs_2, fs_jpim1   ! vector opt. 
    695                   eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk)  
     695                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
    696696               END DO 
    697697            END DO 
     
    720720               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    721721               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    722                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated)  
     722               ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
    723723               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    724724               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 
    725725            END DO 
    726726         END DO 
    727       END DO  
     727      END DO 
    728728 
    729729      ! 
     
    807807            END DO 
    808808         END DO 
     809#if defined key_traldf_c2d || key_traldf_c3d 
    809810         IF( ln_stopack) THEN 
    810811            IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk) 
    811812            IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk) 
    812          ENDIF  
     813         ENDIF 
     814#else 
     815         IF ( ln_stopack ) & 
     816            & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// & 
     817                             'key_traldf_c2d or key_traldf_c3d') 
     818#endif 
    813819      END DO 
    814820      ! 
     
    851857      !!---------------------------------------------------------------------- 
    852858      !!                  ***  ROUTINE zdf_gls_init  *** 
    853       !!                      
    854       !! ** Purpose :   Initialization of the vertical eddy diffivity and  
     859      !! 
     860      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
    855861      !!      viscosity when using a gls turbulent closure scheme 
    856862      !! 
     
    10711077         ! 
    10721078      END SELECT 
    1073      
     1079 
    10741080      !                                !* Set Schmidt number for psi diffusion in the wave breaking case 
    10751081      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 
    10761082      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 
    10771083      IF( ln_sigpsi ) THEN 
    1078          ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf  
     1084         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
    10791085         ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 
    10801086         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 
     
    10841090         rsc_psi0 = rsc_psi 
    10851091      ENDIF 
    1086   
     1092 
    10871093      !                                !* Shear free turbulence parameters 
    10881094      ! 
     
    11301136      rc04  = rc03 * rc0 
    11311137      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1132       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
     1138      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
    11331139      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1134       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1140      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
    11351141      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    1136       rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
     1142      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
    11371143      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
    1138       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
     1144      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
    11391145 
    11401146      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
     
    11511157         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
    11521158      END DO 
    1153       !                               
     1159      ! 
    11541160      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files 
    11551161      ! 
     
    11621168      !!--------------------------------------------------------------------- 
    11631169      !!                   ***  ROUTINE ts_rst  *** 
    1164       !!                      
     1170      !! 
    11651171      !! ** Purpose :   Read or write TKE file (en) in restart file 
    11661172      !! 
    11671173      !! ** Method  :   use of IOM library 
    1168       !!                if the restart does not contain TKE, en is either  
     1174      !!                if the restart does not contain TKE, en is either 
    11691175      !!                set to rn_emin or recomputed (nn_igls/=0) 
    11701176      !!---------------------------------------------------------------------- 
     
    11781184      !!---------------------------------------------------------------------- 
    11791185      ! 
    1180       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     1186      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
    11811187         !                                   ! --------------- 
    11821188         IF( ln_rstart ) THEN                   !* Read the restart file 
     
    11971203               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   ) 
    11981204               IF(nn_timing == 2)  CALL timing_stop('iom_rstget') 
    1199             ELSE                         
     1205            ELSE 
    12001206               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 
    12011207               en  (:,:,:) = rn_emin 
    1202                mxln(:,:,:) = 0.05         
     1208               mxln(:,:,:) = 0.05 
    12031209               avt_k (:,:,:) = avt (:,:,:) 
    12041210               avm_k (:,:,:) = avm (:,:,:) 
     
    12101216            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 
    12111217            en  (:,:,:) = rn_emin 
    1212             mxln(:,:,:) = 0.05        
     1218            mxln(:,:,:) = 0.05 
    12131219         ENDIF 
    12141220         ! 
     
    12171223         IF(lwp) WRITE(numout,*) '---- gls-rst ----' 
    12181224         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
    1219          CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )  
     1225         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
    12201226         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  ) 
    12211227         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  ) 
    1222          CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )  
     1228         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
    12231229         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 
    12241230         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   ) 
     
    12391245      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?' 
    12401246   END SUBROUTINE zdf_gls_init 
    1241     
     1247 
    12421248   SUBROUTINE zdf_gls( kt )          ! Empty routine 
    12431249   IMPLICIT NONE 
    1244       INTEGER, INTENT(in) ::   kt ! ocean time step    
     1250      INTEGER, INTENT(in) ::   kt ! ocean time step 
    12451251      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 
    12461252   END SUBROUTINE zdf_gls 
    1247     
     1253 
    12481254   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine 
    12491255   IMPLICIT NONE 
     
    12561262   !!====================================================================== 
    12571263END MODULE zdfgls 
    1258  
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r11442 r13355  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  zdftke  *** 
    4    !! Ocean physics:  vertical mixing coefficient computed from the tke  
     4   !! Ocean physics:  vertical mixing coefficient computed from the tke 
    55   !!                 turbulent closure parameterization 
    66   !!===================================================================== 
     
    2222   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb 
    2323   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters 
    24    !!             -   !  2008-12  (G. Reffray) stable discretization of the production term  
    25    !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl  
     24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
     25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
    2626   !!                 !                                + cleaning of the parameters + bugs correction 
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     
    5252   USE wrk_nemo       ! work arrays 
    5353   USE timing         ! Timing 
    54    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     54   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    5555#if defined key_agrif 
    5656   USE agrif_opa_interp 
     
    7474   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    7575   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
    76    REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation  
     76   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
    7777   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke 
    7878   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2] 
     
    102102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
    103103   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    104     
     104 
    105105   !! * Substitutions 
    106106#  include "domzgr_substitute.h90" 
     
    118118      !!---------------------------------------------------------------------- 
    119119      ALLOCATE(                                                                    & 
    120          &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &       
     120         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         & 
    121121#if defined key_c1d 
    122122         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          & 
     
    147147      !!         surface: en = max( rn_emin0, rn_ebb * taum ) 
    148148      !!         bottom : en = rn_emin 
    149       !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)  
    150       !! 
    151       !!        The now Turbulent kinetic energy is computed using the following  
     149      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
     150      !! 
     151      !!        The now Turbulent kinetic energy is computed using the following 
    152152      !!      time stepping: implicit for vertical diffusion term, linearized semi 
    153       !!      implicit for kolmogoroff dissipation term, and explicit forward for  
    154       !!      both buoyancy and shear production terms. Therefore a tridiagonal  
     153      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
     154      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
    155155      !!      linear system is solved. Note that buoyancy and shear terms are 
    156156      !!      discretized in a energy conserving form (Bruchard 2002). 
     
    160160      !! 
    161161      !!        The now vertical eddy vicosity and diffusivity coefficients are 
    162       !!      given by:  
     162      !!      given by: 
    163163      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 
    164       !!              avt = max( avmb, pdl * avm                 )   
     164      !!              avt = max( avmb, pdl * avm                 ) 
    165165      !!              eav = max( avmb, avm ) 
    166166      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 
    167       !!      given by an empirical funtion of the localRichardson number if nn_pdl=1  
     167      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
    168168      !! 
    169169      !! ** Action  :   compute en (now turbulent kinetic energy) 
     
    180180      ! 
    181181      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    182          avt (:,:,:) = avt_k (:,:,:)  
    183          avm (:,:,:) = avm_k (:,:,:)  
    184          avmu(:,:,:) = avmu_k(:,:,:)  
    185          avmv(:,:,:) = avmv_k(:,:,:)  
    186       ENDIF 
    187       ! 
     182         avt (:,:,:) = avt_k (:,:,:) 
     183         avm (:,:,:) = avm_k (:,:,:) 
     184         avmu(:,:,:) = avmu_k(:,:,:) 
     185         avmv(:,:,:) = avmv_k(:,:,:) 
     186      ENDIF 
     187      ! 
     188#if defined key_traldf_c2d || key_traldf_c3d 
    188189      IF( ln_stopack) THEN 
    189190         IF( nn_spp_tkelc > 0 ) THEN 
     
    208209         ENDIF 
    209210      ENDIF 
     211#else 
     212      IF ( ln_stopack ) & 
     213         & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// & 
     214                          'key_traldf_c2d or key_traldf_c3d') 
     215#endif 
    210216      ! 
    211217      CALL tke_tke      ! now tke (en) 
     
    213219      CALL tke_avn      ! now avt, avm, avmu, avmv 
    214220      ! 
    215       avt_k (:,:,:) = avt (:,:,:)  
    216       avm_k (:,:,:) = avm (:,:,:)  
    217       avmu_k(:,:,:) = avmu(:,:,:)  
    218       avmv_k(:,:,:) = avmv(:,:,:)  
     221      avt_k (:,:,:) = avt (:,:,:) 
     222      avm_k (:,:,:) = avm (:,:,:) 
     223      avmu_k(:,:,:) = avmu(:,:,:) 
     224      avmv_k(:,:,:) = avmv(:,:,:) 
    219225      ! 
    220226#if defined key_agrif 
    221       ! Update child grid f => parent grid  
     227      ! Update child grid f => parent grid 
    222228      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    223 #endif       
     229#endif 
    224230      IF (  kt == nitend ) THEN 
    225231        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 
     
    271277      CALL wrk_alloc( jpi,jpj, zfact2 ) 
    272278      CALL wrk_alloc( jpi,jpj, zfact3 ) 
    273       CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
     279      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
    274280      ! 
    275281      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation 
    276       zfact1 = -.5_wp * rdt  
     282      zfact1 = -.5_wp * rdt 
    277283      zfact2 = 1.5_wp * rdt * rn_ediss0 
    278284      zfact3 = 0.5_wp       * rn_ediss0 
     
    294300         END DO 
    295301      END DO 
    296        
     302 
    297303!!bfr   - start commented area 
    298304      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    333339         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    334340         DO jk = jpkm1, 2, -1 
    335             DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us  
     341            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
    336342               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1) 
    337343                  zus  = zcof * taum(ji,jj) 
     
    341347         END DO 
    342348         !                               ! finite LC depth 
    343          DO jj = 1, jpj  
     349         DO jj = 1, jpj 
    344350            DO ji = 1, jpi 
    345351               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 
     
    378384            DO ji = 1, jpi 
    379385               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    380                   &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
     386                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
    381387                  &                            / (  fse3uw_n(ji,jj,jk)               & 
    382388                  &                              *  fse3uw_b(ji,jj,jk)  ) 
     
    407413                  !                                                           ! shear prod. at w-point weightened by mask 
    408414               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    409                   &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
     415                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 
    410416                  ! 
    411417               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
     
    464470      END DO 
    465471 
    466       !                                 ! Save TKE prior to nn_etau addition   
    467       e_niw(:,:,:) = en(:,:,:)   
    468       !   
     472      !                                 ! Save TKE prior to nn_etau addition 
     473      e_niw(:,:,:) = en(:,:,:) 
     474      ! 
    469475      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    470476      !                            !  TKE due to surface and internal wave breaking 
     
    508514                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    509515                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    510                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    511                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
     516                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress 
     517                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
    512518                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    513519                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     
    517523         END DO 
    518524      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up) 
    519          IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau  
     525         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
    520526            DO jj = 2, jpjm1 
    521527               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    535541      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    536542      ! 
    537       DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term   
    538          DO jj = 2, jpjm1   
    539             DO ji = fs_2, fs_jpim1   ! vector opt.   
    540                e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)   
    541             END DO   
    542          END DO   
    543       END DO   
    544       !   
    545       CALL lbc_lnk( e_niw, 'W', 1. )   
     543      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term 
     544         DO jj = 2, jpjm1 
     545            DO ji = fs_2, fs_jpim1   ! vector opt. 
     546               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 
     547            END DO 
     548         END DO 
     549      END DO 
     550      ! 
     551      CALL lbc_lnk( e_niw, 'W', 1. ) 
    546552      ! 
    547553      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    548       CALL wrk_dealloc( jpi,jpj, zhlc )  
    549       CALL wrk_dealloc( jpi,jpj, zbbrau )  
    550       CALL wrk_dealloc( jpi,jpj, zfact2 )  
    551       CALL wrk_dealloc( jpi,jpj, zfact3 )  
    552       CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
     554      CALL wrk_dealloc( jpi,jpj, zhlc ) 
     555      CALL wrk_dealloc( jpi,jpj, zbbrau ) 
     556      CALL wrk_dealloc( jpi,jpj, zfact2 ) 
     557      CALL wrk_dealloc( jpi,jpj, zfact3 ) 
     558      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
    553559      ! 
    554560      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
     
    563569      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity 
    564570      !! 
    565       !! ** Method  :   At this stage, en, the now TKE, is known (computed in  
    566       !!              the tke_tke routine). First, the now mixing lenth is  
     571      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
     572      !!              the tke_tke routine). First, the now mixing lenth is 
    567573      !!      computed from en and the strafification (N^2), then the mixings 
    568574      !!      coefficients are computed. 
    569575      !!              - Mixing length : a first evaluation of the mixing lengh 
    570576      !!      scales is: 
    571       !!                      mxl = sqrt(2*en) / N   
     577      !!                      mxl = sqrt(2*en) / N 
    572578      !!      where N is the brunt-vaisala frequency, with a minimum value set 
    573579      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean. 
    574       !!        The mixing and dissipative length scale are bound as follow :  
     580      !!        The mixing and dissipative length scale are bound as follow : 
    575581      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom. 
    576582      !!                        zmxld = zmxlm = mxl 
    577583      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 
    578       !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is  
     584      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
    579585      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 
    580586      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings 
    581       !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to  
    582       !!                    the surface to obtain ldown. the resulting length  
     587      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
     588      !!                    the surface to obtain ldown. the resulting length 
    583589      !!                    scales are: 
    584       !!                        zmxld = sqrt( lup * ldown )  
     590      !!                        zmxld = sqrt( lup * ldown ) 
    585591      !!                        zmxlm = min ( lup , ldown ) 
    586592      !!              - Vertical eddy viscosity and diffusivity: 
    587593      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 
    588       !!                      avt = max( avmb, pdlr * avm )   
     594      !!                      avt = max( avmb, pdlr * avm ) 
    589595      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 
    590596      !! 
     
    601607      IF( nn_timing == 1 )  CALL timing_start('tke_avn') 
    602608 
    603       CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
     609      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
    604610 
    605611      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    610616      ! 
    611617      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
    612       zmxlm(:,:,:)  = rmxl_min     
     618      zmxlm(:,:,:)  = rmxl_min 
    613619      zmxld(:,:,:)  = rmxl_min 
    614620      ! 
     
    620626            END DO 
    621627         END DO 
    622       ELSE  
     628      ELSE 
    623629         zmxlm(:,:,1) = rn_mxl0 
    624630      ENDIF 
     
    638644      !                     !* Physical limits for the mixing length 
    639645      ! 
    640       zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value  
     646      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
    641647      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    642648      ! 
     
    738744            END DO 
    739745         END DO 
     746#if defined key_traldf_c2d || key_traldf_c3d 
    740747         IF( ln_stopack) THEN 
    741748            IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 
    742749            IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 
    743750         ENDIF 
     751#else 
     752         IF ( ln_stopack  ) & 
     753            & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// & 
     754                             'key_traldf_c2d or key_traldf_c3d') 
     755#endif 
    744756      END DO 
    745757      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     
    787799      ENDIF 
    788800      ! 
    789       CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
     801      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
    790802      ! 
    791803      IF( nn_timing == 1 )  CALL timing_stop('tke_avn') 
     
    797809      !!---------------------------------------------------------------------- 
    798810      !!                  ***  ROUTINE zdf_tke_init  *** 
    799       !!                      
    800       !! ** Purpose :   Initialization of the vertical eddy diffivity and  
     811      !! 
     812      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
    801813      !!              viscosity when using a tke turbulent closure scheme 
    802814      !! 
     
    814826         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   & 
    815827         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   & 
    816          &                 nn_etau , nn_htau  , rn_efr , rn_c    
     828         &                 nn_etau , nn_htau  , rn_efr , rn_c 
    817829      !!---------------------------------------------------------------------- 
    818830 
     
    884896      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb 
    885897      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr 
    886        
     898 
    887899      IF( nn_etau == 2  ) THEN 
    888900          ierr = zdf_mxl_alloc() 
     
    896908 
    897909      !                               !* depth of penetration of surface tke 
    898       IF( nn_etau /= 0 ) THEN       
     910      IF( nn_etau /= 0 ) THEN 
    899911         htau(:,:) = 0._wp 
    900912         SELECT CASE( nn_htau )             ! Choice of the depth of penetration 
     
    902914            htau(:,:) = 10._wp 
    903915         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees 
    904             htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )             
     916            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   ) 
    905917         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer 
    906918            rhtau = -1._wp / LOG( 1._wp - rn_c ) 
     
    945957      END DO 
    946958      dissl(:,:,:) = 1.e-12_wp 
    947       !                               
     959      ! 
    948960      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files 
    949961      ! 
     
    954966     !!--------------------------------------------------------------------- 
    955967     !!                   ***  ROUTINE tke_rst  *** 
    956      !!                      
     968     !! 
    957969     !! ** Purpose :   Read or write TKE file (en) in restart file 
    958970     !! 
    959971     !! ** Method  :   use of IOM library 
    960      !!                if the restart does not contain TKE, en is either  
    961      !!                set to rn_emin or recomputed  
     972     !!                if the restart does not contain TKE, en is either 
     973     !!                set to rn_emin or recomputed 
    962974     !!---------------------------------------------------------------------- 
    963975     INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     
    968980     !!---------------------------------------------------------------------- 
    969981     ! 
    970      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     982     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
    971983        !                                   ! --------------- 
    972984        IF( ln_rstart ) THEN                   !* Read the restart file 
     
    10061018        avmu_k(:,:,:) = avmu(:,:,:) 
    10071019        avmv_k(:,:,:) = avmv(:,:,:) 
    1008          
     1020 
    10091021     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    10101022        !                                   ! ------------------- 
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/step.F90

    r11442 r13355  
    110110      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    111111      IF( lk_bdy     )  THEN 
    112          IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib  
     112         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
    113113                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    114114      ENDIF 
     
    118118         CALL lbc_lnk( tsb(:,:,:,tind), 'T', 1. ) 
    119119      END DO 
    120        
     120 
    121121      IF( ln_stopack )   CALL stopack_pert( kstp ) 
    122122                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     
    153153      ENDIF 
    154154      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
     155#if defined key_traldf_c2d || key_traldf_c3d 
    155156         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) THEN 
    156157              rn_avt_rnf0 = rn_avt_rnf 
    157158              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf ) 
    158159         ENDIF 
     160#else 
     161         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) & 
     162            & CALL ctl_stop( 'stp: parameter perturbation will only work with '// & 
     163                             'key_traldf_c2d or key_traldf_c3d') 
     164#endif 
    159165         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
    160166      ENDIF 
     
    204210      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    205211                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur) 
    206       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    207                          CALL wzv           ( kstp )  ! now cross-level velocity  
    208  
    209       IF( lk_dynspg_ts ) THEN  
     212      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
     213                         CALL wzv           ( kstp )  ! now cross-level velocity 
     214 
     215      IF( lk_dynspg_ts ) THEN 
    210216          ! In case the time splitting case, update almost all momentum trends here: 
    211217          ! Note that the computation of vertical velocity above, hence "after" sea level 
     
    244250                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 
    245251          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    246                                   CALL wzv           ( kstp )  ! now cross-level velocity  
     252                                  CALL wzv           ( kstp )  ! now cross-level velocity 
    247253      ENDIF 
    248254 
     
    320326               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    321327               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    322          IF( ln_zps .AND.       ln_isfcav)                                   &  
     328         IF( ln_zps .AND.       ln_isfcav)                                   & 
    323329               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    324330               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     
    334340                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    335341                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    336                                 IF( ln_zps .AND.       ln_isfcav)                                   &  
     342                                IF( ln_zps .AND.       ln_isfcav)                                   & 
    337343                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    338344                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     
    349355                               ua(:,:,:) = ua_sv(:,:,:) 
    350356                               va(:,:,:) = va_sv(:,:,:) 
    351                                                              ! Revert now divergence and rotational to previously computed ones  
     357                                                             ! Revert now divergence and rotational to previously computed ones 
    352358                                                             !(needed because of the time swap in div_cur, at the beginning of each time step) 
    353359                               hdivn(:,:,:) = hdivb(:,:,:) 
    354                                rotn(:,:,:)  = rotb(:,:,:)  
     360                               rotn(:,:,:)  = rotb(:,:,:) 
    355361 
    356362                               CALL dyn_bfr( kstp )         ! bottom friction 
     
    398404      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    399405      ! AGRIF 
    400       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    401                                CALL Agrif_Integrate_ChildGrids( stp )   
     406      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     407                               CALL Agrif_Integrate_ChildGrids( stp ) 
    402408 
    403409      IF ( Agrif_NbStepint().EQ.0 ) THEN 
     
    431437      ! 
    432438#if defined key_iomput 
    433       IF( kstp == nitend .OR. indic < 0 ) THEN  
     439      IF( kstp == nitend .OR. indic < 0 ) THEN 
    434440                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    435          IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !  
     441         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
    436442      ENDIF 
    437443#endif 
    438444      ! 
    439445      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset 
    440       !      
     446      ! 
    441447      ! 
    442448   END SUBROUTINE stp 
Note: See TracChangeset for help on using the changeset viewer.