New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13088 for branches/UKMO – NEMO

Changeset 13088 for branches/UKMO


Ignore:
Timestamp:
2020-06-10T13:13:39+02:00 (4 years ago)
Author:
jwhile
Message:

Bug fixes for 1D running

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

Legend:

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

    r10302 r13088  
    3333   USE oce, ONLY:           & ! active tracer variables 
    3434      & tsn 
    35    USE zdfmxl, ONLY :       & ! mixed layer depth  
     35   USE zdfmxl, ONLY :       & ! mixed layer depth 
    3636#if defined key_karaml 
    3737      & hmld_kara,          & 
    3838      & ln_kara,            & 
    39 #endif    
    40       & hmld,               &  
     39#endif 
     40      & hmld,               & 
    4141      & hmlp,               & 
    42       & hmlpt  
     42      & hmlpt 
    4343   USE asmpar, ONLY:        & ! assimilation parameters 
    4444      & c_asmbkg,           & 
     
    8989 
    9090   IMPLICIT NONE 
    91    PRIVATE                    
     91   PRIVATE 
    9292 
    9393   PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90 
     
    290290 
    291291      ! Allocate and read increments 
    292        
     292 
    293293      IF ( ln_slchltotinc ) THEN 
    294294         ALLOCATE( slchltot_bkginc(jpi,jpj) ) 
    295295         CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc ) 
    296296      ENDIF 
    297        
     297 
    298298      IF ( ln_slchldiainc ) THEN 
    299299         ALLOCATE( slchldia_bkginc(jpi,jpj) ) 
    300300         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc ) 
    301301      ENDIF 
    302        
     302 
    303303      IF ( ln_slchlnoninc ) THEN 
    304304         ALLOCATE( slchlnon_bkginc(jpi,jpj) ) 
    305305         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc ) 
    306306      ENDIF 
    307        
     307 
    308308      IF ( ln_schltotinc ) THEN 
    309309         ALLOCATE( schltot_bkginc(jpi,jpj) ) 
    310310         CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc ) 
    311311      ENDIF 
    312        
     312 
    313313      IF ( ln_slphytotinc ) THEN 
    314314         ALLOCATE( slphytot_bkginc(jpi,jpj) ) 
    315315         CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc ) 
    316316      ENDIF 
    317        
     317 
    318318      IF ( ln_slphydiainc ) THEN 
    319319         ALLOCATE( slphydia_bkginc(jpi,jpj) ) 
    320320         CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc ) 
    321321      ENDIF 
    322        
     322 
    323323      IF ( ln_slphynoninc ) THEN 
    324324         ALLOCATE( slphynon_bkginc(jpi,jpj) ) 
     
    335335         CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc ) 
    336336      ENDIF 
    337        
     337 
    338338      IF ( ln_plchltotinc ) THEN 
    339339         ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) ) 
    340340         CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc ) 
    341341      ENDIF 
    342        
     342 
    343343      IF ( ln_pchltotinc ) THEN 
    344344         ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) ) 
    345345         CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc ) 
    346346      ENDIF 
    347        
     347 
    348348      IF ( ln_pno3inc ) THEN 
    349349         ALLOCATE( pno3_bkginc(jpi,jpj,jpk) ) 
    350350         CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc ) 
    351351      ENDIF 
    352        
     352 
    353353      IF ( ln_psi4inc ) THEN 
    354354         ALLOCATE( psi4_bkginc(jpi,jpj,jpk) ) 
    355355         CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc ) 
    356356      ENDIF 
    357        
     357 
    358358      IF ( ln_pdicinc ) THEN 
    359359         ALLOCATE( pdic_bkginc(jpi,jpj,jpk) ) 
    360360         CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc ) 
    361361      ENDIF 
    362        
     362 
    363363      IF ( ln_palkinc ) THEN 
    364364         ALLOCATE( palk_bkginc(jpi,jpj,jpk) ) 
    365365         CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc ) 
    366366      ENDIF 
    367        
     367 
    368368      IF ( ln_pphinc ) THEN 
    369369         ALLOCATE( pph_bkginc(jpi,jpj,jpk) ) 
    370370         CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc ) 
    371371      ENDIF 
    372        
     372 
    373373      IF ( ln_po2inc ) THEN 
    374374         ALLOCATE( po2_bkginc(jpi,jpj,jpk) ) 
     
    377377 
    378378      ! Allocate balancing increments 
    379        
     379 
    380380      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
    381381         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
     
    388388#endif 
    389389      ENDIF 
    390        
     390 
    391391      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN 
    392392#if defined key_top 
     
    443443      ! Initialise 
    444444      p_incs(:,:) = 0.0 
    445        
     445 
    446446      ! read from file 
    447447      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 ) 
    448        
     448 
    449449      ! Apply the masks 
    450450      p_incs(:,:) = p_incs(:,:) * tmask(:,:,1) 
    451        
     451 
    452452      ! Set missing increments to 0.0 rather than 1e+20 
    453453      ! to allow for differences in masks 
     
    481481      ! Initialise 
    482482      p_incs(:,:,:) = 0.0 
    483        
     483 
    484484      ! read from file 
    485485      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 ) 
    486        
     486 
    487487      ! Apply the masks 
    488488      p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:) 
    489        
     489 
    490490      ! Set missing increments to 0.0 rather than 1e+20 
    491491      ! to allow for differences in masks 
     
    538538         cchl_p_bkg(:,:,:) = 0.0 
    539539#endif 
    540           
     540 
    541541         !-------------------------------------------------------------------- 
    542542         ! Read background variables for phytoplankton assimilation 
     
    558558         CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 
    559559#endif 
    560           
     560 
    561561         IF ( ln_phytobal ) THEN 
    562562 
     
    602602 
    603603         CALL iom_close( inum ) 
    604           
     604 
    605605         DO jt = 1, jptra 
    606606            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    607607         END DO 
    608        
     608 
    609609      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN 
    610610 
     
    615615 
    616616         CALL iom_open( c_asmbkg, inum ) 
    617           
     617 
    618618#if defined key_hadocc 
    619619         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 
     
    626626 
    627627         CALL iom_close( inum ) 
    628           
     628 
    629629         DO jt = 1, jptra 
    630630            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 
    631631         END DO 
    632632         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 
    633        
     633 
    634634      ENDIF 
    635635#else 
     
    655655      !! 
    656656      !! ** Action  : 
    657       !!                    
     657      !! 
    658658      !! References : 
    659659      !! 
     
    669669      REAL(wp)          :: zdate         ! Date 
    670670      !!------------------------------------------------------------------------ 
    671       
     671 
    672672      ! Set things up 
    673673      zdate = REAL( ndastp ) 
     
    680680            &                    TRIM( c_asmbal ) // ' at timestep = ', kt 
    681681 
    682          ! Define the output file        
     682         ! Define the output file 
    683683         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib) 
    684684 
     
    767767            &                    TRIM( c_asmbal ) // ' at timestep = ', kt 
    768768      ENDIF 
    769                                   
     769 
    770770   END SUBROUTINE asm_bgc_bal_wri 
    771771 
     
    847847      !! ** Action  :   return non-log increments 
    848848      !! 
    849       !! References :    
     849      !! References : 
    850850      !!------------------------------------------------------------------------ 
    851851      !! 
     
    879879      !!------------------------------------------------------------------------ 
    880880      !!                    ***  ROUTINE phyto2d_asm_inc  *** 
    881       !!           
     881      !! 
    882882      !! ** Purpose : Apply the chlorophyll assimilation increments. 
    883883      !! 
     
    886886      !!              Direct initialization or Incremental Analysis Updating. 
    887887      !! 
    888       !! ** Action  :  
     888      !! ** Action  : 
    889889      !!------------------------------------------------------------------------ 
    890890      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    914914#endif 
    915915      !!------------------------------------------------------------------------ 
    916        
     916 
    917917      IF ( kt <= nit000 ) THEN 
    918918 
     
    920920         ! Remember that two sets of non-log increments should not be 
    921921         ! expected to be in the same ratio as their log equivalents 
    922           
     922 
    923923         ! Total chlorophyll 
    924924         IF ( ln_slchltotinc ) THEN 
     
    10741074 
    10751075            IF(lwp) THEN 
    1076                WRITE(numout,*)  
     1076               WRITE(numout,*) 
    10771077               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', & 
    10781078                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    11051105         ENDIF 
    11061106 
    1107       ELSEIF ( ll_asmdin ) THEN  
     1107      ELSEIF ( ll_asmdin ) THEN 
    11081108 
    11091109         !-------------------------------------------------------------------- 
    11101110         ! Direct Initialization 
    11111111         !-------------------------------------------------------------------- 
    1112           
     1112 
    11131113         IF ( kt == nitdin_r ) THEN 
    11141114 
     
    11341134            END WHERE 
    11351135#endif 
    1136   
     1136 
    11371137            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    11381138            ! which is called at end of model run 
     
    11501150      !!------------------------------------------------------------------------ 
    11511151      !!                    ***  ROUTINE phyto3d_asm_inc  *** 
    1152       !!           
     1152      !! 
    11531153      !! ** Purpose : Apply the profile chlorophyll assimilation increments. 
    11541154      !! 
     
    11561156      !!              Direct initialization or Incremental Analysis Updating. 
    11571157      !! 
    1158       !! ** Action  :  
     1158      !! ** Action  : 
    11591159      !!------------------------------------------------------------------------ 
    11601160      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    12611261 
    12621262            IF(lwp) THEN 
    1263                WRITE(numout,*)  
     1263               WRITE(numout,*) 
    12641264               WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', & 
    12651265                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    12921292         ENDIF 
    12931293 
    1294       ELSEIF ( ll_asmdin ) THEN  
     1294      ELSEIF ( ll_asmdin ) THEN 
    12951295 
    12961296         !-------------------------------------------------------------------- 
    12971297         ! Direct Initialization 
    12981298         !-------------------------------------------------------------------- 
    1299           
     1299 
    13001300         IF ( kt == nitdin_r ) THEN 
    13011301 
     
    13211321            END WHERE 
    13221322#endif 
    1323   
     1323 
    13241324            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    13251325            ! which is called at end of model run 
     
    13381338      !!------------------------------------------------------------------------ 
    13391339      !!                    ***  ROUTINE pco2_asm_inc  *** 
    1340       !!           
     1340      !! 
    13411341      !! ** Purpose : Apply the pco2/fco2 assimilation increments. 
    13421342      !! 
     
    13451345      !!              Direct initialization or Incremental Analysis Updating. 
    13461346      !! 
    1347       !! ** Action  :  
     1347      !! ** Action  : 
    13481348      !!------------------------------------------------------------------------ 
    13491349      INTEGER, INTENT(IN)                   :: kt           ! Current time step 
     
    14951495               jkmax = jpk-1 
    14961496               DO jk = jpk-1, 1, -1 
     1497#if defined key_vvl 
    14971498                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
    14981499                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
     
    15001501                     jkmax = jk 
    15011502                  ENDIF 
     1503#else 
     1504                  IF ( ( zmld(ji,jj) >  gdepw_0(ji,jj,jk)   ) .AND. & 
     1505                     & ( zmld(ji,jj) <= gdepw_0(ji,jj,jk+1) ) ) THEN 
     1506                     zmld(ji,jj) = gdepw_0(ji,jj,jk+1) 
     1507                     jkmax = jk 
     1508                  ENDIF 
     1509#endif 
    15021510               END DO 
    15031511               ! 
     
    15261534 
    15271535            IF(lwp) THEN 
    1528                WRITE(numout,*)  
     1536               WRITE(numout,*) 
    15291537               IF ( ln_spco2inc ) THEN 
    15301538                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 
     
    15631571         ENDIF 
    15641572 
    1565       ELSEIF ( ll_asmdin ) THEN  
     1573      ELSEIF ( ll_asmdin ) THEN 
    15661574 
    15671575         !-------------------------------------------------------------------- 
    15681576         ! Direct Initialization 
    15691577         !-------------------------------------------------------------------- 
    1570           
     1578 
    15711579         IF ( kt == nitdin_r ) THEN 
    15721580 
     
    15921600            END WHERE 
    15931601#endif 
    1594   
     1602 
    15951603            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    15961604            ! which is called at end of model run 
     
    16091617      !!------------------------------------------------------------------------ 
    16101618      !!                    ***  ROUTINE ph_asm_inc  *** 
    1611       !!           
     1619      !! 
    16121620      !! ** Purpose : Apply the pH assimilation increments. 
    16131621      !! 
     
    16161624      !!              Direct initialization or Incremental Analysis Updating. 
    16171625      !! 
    1618       !! ** Action  :  
     1626      !! ** Action  : 
    16191627      !!------------------------------------------------------------------------ 
    16201628      INTEGER,                          INTENT(IN) :: kt        ! Current time step 
     
    16261634      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments 
    16271635      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments 
    1628        
     1636 
    16291637      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation 
    16301638      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation 
     
    16791687            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) 
    16801688         ENDIF 
    1681           
     1689 
    16821690         ! Account for pCO2 balancing if required 
    16831691         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN 
     
    16851693            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) 
    16861694         ENDIF 
    1687           
     1695 
    16881696         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk 
    16891697         ! This requires three calls to the MOCSY carbonate package 
     
    17501758                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic 
    17511759                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk 
    1752                       
     1760 
    17531761                  ENDIF 
    1754                    
     1762 
    17551763               END DO 
    17561764            END DO 
     
    17581766 
    17591767      ENDIF 
    1760        
     1768 
    17611769      IF ( ll_asmiau ) THEN 
    17621770 
     
    17721780 
    17731781            IF(lwp) THEN 
    1774                WRITE(numout,*)  
     1782               WRITE(numout,*) 
    17751783               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & 
    17761784                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    17941802         ENDIF 
    17951803 
    1796       ELSEIF ( ll_asmdin ) THEN  
     1804      ELSEIF ( ll_asmdin ) THEN 
    17971805 
    17981806         !-------------------------------------------------------------------- 
    17991807         ! Direct Initialization 
    18001808         !-------------------------------------------------------------------- 
    1801           
     1809 
    18021810         IF ( kt == nitdin_r ) THEN 
    18031811 
     
    18141822               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
    18151823            END WHERE 
    1816   
     1824 
    18171825            ! Do not deallocate arrays - needed by asm_bgc_bal_wri 
    18181826            ! which is called at end of model run 
     
    18201828         ! 
    18211829      ENDIF 
    1822 #endif       
     1830#endif 
    18231831      ! 
    18241832   END SUBROUTINE ph_asm_inc 
     
    18311839      !!---------------------------------------------------------------------- 
    18321840      !!                    ***  ROUTINE dyn_asm_inc  *** 
    1833       !!           
     1841      !! 
    18341842      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments. 
    18351843      !! 
    18361844      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
    18371845      !! 
    1838       !! ** Action  :  
     1846      !! ** Action  : 
    18391847      !!---------------------------------------------------------------------- 
    18401848      INTEGER,  INTENT(IN) :: kt        ! Current time step 
     
    19831991 
    19841992            IF(lwp) THEN 
    1985                WRITE(numout,*)  
     1993               WRITE(numout,*) 
    19861994               WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', & 
    19871995                  &  kt,' with IAU weight = ', pwgtiau(it) 
     
    20712079#endif 
    20722080            ENDIF 
    2073             
     2081 
    20742082            IF ( kt == nitiaufin_r ) THEN 
    20752083               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) 
     
    20822090         ENDIF 
    20832091 
    2084       ELSEIF ( ll_asmdin ) THEN  
     2092      ELSEIF ( ll_asmdin ) THEN 
    20852093 
    20862094         !-------------------------------------------------------------------- 
    20872095         ! Direct Initialization 
    20882096         !-------------------------------------------------------------------- 
    2089           
     2097 
    20902098         IF ( kt == nitdin_r ) THEN 
    20912099 
     
    21792187#endif 
    21802188            ENDIF 
    2181   
     2189 
    21822190            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc ) 
    21832191            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc ) 
  • branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r8400 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90

    r12102 r13088  
    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, 
     
    624624      !!---------------------------------------------------------------------- 
    625625   INTEGER, INTENT( in ) :: kt 
     626#if defined key_traldf_c3d 
     627   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff 
     628   REAL(wp), POINTER, DIMENSION(:,:,:) ::   gauss 
     629#elif defined key_traldf_c2d 
    626630   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff 
     631   REAL(wp), POINTER, DIMENSION(:,:) ::   gauss 
     632#elif defined key_traldf_c1d 
     633   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff 
     634   REAL(wp), POINTER, DIMENSION(:) ::   gauss 
     635#else 
     636   REAL(wp), INTENT( inout ) :: coeff 
     637   REAL(wp), POINTER ::   gauss 
     638#endif 
    627639   INTEGER, INTENT( in ) ::  nn_type 
    628640   REAL(wp), INTENT( in ) :: rn_sd 
     
    634646   INTEGER :: jklev 
    635647 
     648#if defined key_traldf_c2d || key_traldf_c3d 
    636649   CALL wrk_alloc(jpi,jpj,gauss) 
    637650 
     
    658671       gauss = gauss * rn_sd 
    659672       coeff = coeff * ( 1._wp + gauss ) 
     673#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d 
    660674       WHERE( coeff > 1._wp ) coeff = 1._wp 
    661675       WHERE( coeff < 0._wp ) coeff = 0._wp 
     676#else 
     677       IF( coeff > 1._wp ) coeff = 1._wp 
     678       IF( coeff < 0._wp ) coeff = 0._wp 
     679#endif 
    662680   ELSEIF ( nn_type  == 5 ) THEN 
    663681       zsd = rn_sd 
     
    665683       gauss = gauss * zsd + xme 
    666684       coeff = exp(gauss) * coeff 
     685#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d 
    667686       WHERE( coeff > 1._wp ) coeff = 1._wp 
    668687       WHERE( coeff < 0._wp ) coeff = 0._wp 
     688#else 
     689       IF( coeff > 1._wp ) coeff = 1._wp 
     690       IF( coeff < 0._wp ) coeff = 0._wp 
     691#endif 
    669692   ELSEIF ( nn_type == 6 ) THEN 
    670693       zsd = rn_sd 
     
    672695       gauss = gauss * zsd + xme 
    673696       coeff = exp(gauss) * coeff 
     697#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d 
    674698       WHERE( coeff > 1._wp ) coeff = 1._wp 
    675699       WHERE( coeff < 0._wp ) coeff = 0._wp 
     700#else 
     701       IF( coeff > 1._wp ) coeff = 1._wp 
     702       IF( coeff < 0._wp ) coeff = 0._wp 
     703#endif 
    676704   ELSE 
    677705       CALL ctl_stop( 'spp dqdt wrong option') 
     
    687715       jklev = klev 
    688716     ELSE 
    689        jklev = 0  
     717       jklev = 0 
    690718     ENDIF 
    691719     CALL spp_stats(kt,kspp,jklev,coeff) 
     
    694722   CALL wrk_dealloc(jpi,jpj,gauss) 
    695723 
     724#else 
     725   CALL ctl_stop( 'spp_gen: parameter perturbation will only work with '// & 
     726                  'key_traldf_c2d or key_traldf_c3d') 
     727#endif 
     728 
     729 
    696730   END SUBROUTINE 
    697731 
     
    704738   IMPLICIT NONE 
    705739   INTEGER, INTENT(IN) :: mt,kp,kl 
    706    REAL(wp), INTENT(IN) :: rcf(jpi,jpj) 
     740#if defined key_traldf_c3d 
     741   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: rcf 
     742#elif defined key_traldf_c2d 
     743   REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: rcf 
     744#elif defined key_traldf_c1d 
     745   REAL(wp), INTENT( inout ), DIMENSION(jpk) :: rcf 
     746#else 
     747   REAL(wp), INTENT( inout ) :: rcf 
     748#endif 
    707749   REAL(wp) :: mi,ma 
    708750   CHARACTER(LEN=16) :: cstr = '                ' 
    709    SELECT CASE ( kp )  
    710      CASE( jk_spp_alb   )  
     751   SELECT CASE ( kp ) 
     752     CASE( jk_spp_alb   ) 
    711753       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   )  
     754     CASE( jk_spp_rhg   ) 
     755       cstr = 'ICE RHEOLOGY' 
     756     CASE( jk_spp_relw  ) 
     757       cstr = 'RELATIVE WND' 
     758     CASE( jk_spp_dqdt  ) 
     759       cstr = 'SST RELAXAT.' 
     760     CASE( jk_spp_deds  ) 
     761       cstr = 'SSS RELAXAT.' 
     762     CASE( jk_spp_arnf  ) 
     763       cstr = 'RIVER MIXING' 
     764     CASE( jk_spp_geot  ) 
     765       cstr = 'GEOTHERM.FLX' 
     766     CASE( jk_spp_qsi0  ) 
     767       cstr = 'SOLAR EXTIN.' 
     768     CASE( jk_spp_bfr   ) 
     769       cstr = 'BOTTOM FRICT' 
     770     CASE( jk_spp_aevd  ) 
     771       cstr = 'EDDY VISCDIF' 
     772     CASE( jk_spp_avt   ) 
    731773       cstr = 'VERT. DIFFUS' 
    732      CASE( jk_spp_avm   )  
     774     CASE( jk_spp_avm   ) 
    733775       cstr = 'VERT. VISCOS' 
    734      CASE( jk_spp_tkelc )  
     776     CASE( jk_spp_tkelc ) 
    735777       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 )  
     778     CASE( jk_spp_tkedf ) 
     779       cstr = 'TKE RN_EDIFF' 
     780     CASE( jk_spp_tkeds ) 
     781       cstr = 'TKE RN_EDISS' 
     782     CASE( jk_spp_tkebb ) 
    741783       cstr = 'TKE RN_EBB  ' 
    742      CASE( jk_spp_tkefr )  
     784     CASE( jk_spp_tkefr ) 
    743785       cstr = 'TKE RN_EFR  ' 
    744      CASE( jk_spp_ahtu  )  
     786     CASE( jk_spp_ahtu  ) 
    745787       cstr = 'TRALDF AHTU ' 
    746      CASE( jk_spp_ahtv  )  
     788     CASE( jk_spp_ahtv  ) 
    747789       cstr = 'TRALDF AHTV ' 
    748      CASE( jk_spp_ahtw  )  
     790     CASE( jk_spp_ahtw  ) 
    749791       cstr = 'TRALDF AHTW ' 
    750      CASE( jk_spp_ahtt  )  
     792     CASE( jk_spp_ahtt  ) 
    751793       cstr = 'TRALDF AHTT ' 
    752794     CASE( jk_spp_ahubbl ) 
     
    765807       CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics') 
    766808   END SELECT 
     809#ifdef key_traldf_c3d || key_traldf_c2d || key_traldf_c1d 
    767810   mi = MINVAL(rcf) 
    768811   ma = MAXVAL(rcf) 
     812#else 
     813   mi = rcf 
     814   ma = rcf 
     815#endif 
    769816   IF(lk_mpp) CALL mpp_min(mi) 
    770817   IF(lk_mpp) CALL mpp_max(ma) 
     
    795842 
    796843   DO jp=1,jk_spp 
    797      SELECT CASE ( jp )  
    798        CASE( jk_spp_alb   )  
     844     SELECT CASE ( jp ) 
     845       CASE( jk_spp_alb   ) 
    799846         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   )  
     847       CASE( jk_spp_rhg   ) 
     848         cstr = 'ICE RHEOLOGY' 
     849       CASE( jk_spp_relw  ) 
     850         cstr = 'RELATIVE WND' 
     851       CASE( jk_spp_dqdt  ) 
     852         cstr = 'SST RELAXAT.' 
     853       CASE( jk_spp_deds  ) 
     854         cstr = 'SSS RELAXAT.' 
     855       CASE( jk_spp_arnf  ) 
     856         cstr = 'RIVER MIXING' 
     857       CASE( jk_spp_geot  ) 
     858         cstr = 'GEOTHERM.FLX' 
     859       CASE( jk_spp_qsi0  ) 
     860         cstr = 'SOLAR EXTIN.' 
     861       CASE( jk_spp_bfr   ) 
     862         cstr = 'BOTTOM FRICT' 
     863       CASE( jk_spp_aevd  ) 
     864         cstr = 'EDDY VISCDIF' 
     865       CASE( jk_spp_avt   ) 
    819866         cstr = 'VERT. DIFFUS' 
    820        CASE( jk_spp_avm   )  
     867       CASE( jk_spp_avm   ) 
    821868         cstr = 'VERT. VISCOS' 
    822        CASE( jk_spp_tkelc )  
     869       CASE( jk_spp_tkelc ) 
    823870         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 )  
     871       CASE( jk_spp_tkedf ) 
     872         cstr = 'TKE RN_EDIFF' 
     873       CASE( jk_spp_tkeds ) 
     874         cstr = 'TKE RN_EDISS' 
     875       CASE( jk_spp_tkebb ) 
    829876         cstr = 'TKE RN_EBB  ' 
    830        CASE( jk_spp_tkefr )  
     877       CASE( jk_spp_tkefr ) 
    831878         cstr = 'TKE RN_EFR  ' 
    832        CASE( jk_spp_ahtu  )  
     879       CASE( jk_spp_ahtu  ) 
    833880         cstr = 'TRALDF AHTU ' 
    834        CASE( jk_spp_ahtv  )  
     881       CASE( jk_spp_ahtv  ) 
    835882         cstr = 'TRALDF AHTV ' 
    836        CASE( jk_spp_ahtw  )  
     883       CASE( jk_spp_ahtw  ) 
    837884         cstr = 'TRALDF AHTW ' 
    838        CASE( jk_spp_ahtt  )  
     885       CASE( jk_spp_ahtt  ) 
    839886         cstr = 'TRALDF AHTT ' 
    840887       CASE( jk_spp_ahubbl ) 
     
    11951242      ! Note: sshn should be staggered before being used. 
    11961243      SELECT CASE ( cd_type ) 
    1197                CASE ( 'T' )   
     1244               CASE ( 'T' ) 
    11981245                jk=1 
    11991246                zv = SUM( tmask_i(:,:)*tmask(:,:,jk)*e1t(:,:)*e2t(:,:)*sshn(:,:)*zts(:,:,1) ) 
     
    12851332      ! Random noise on 2d field 
    12861333      IF ( istep == 1 ) THEN 
    1287          CALL sppt_rand2d( g2d )  
     1334         CALL sppt_rand2d( g2d ) 
    12881335         CALL lbc_lnk( g2d , 'T', 1._wp) 
    12891336         g2d_save = g2d 
     
    12971344         g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev 
    12981345      ENDIF 
    1299     
     1346 
    13001347      ! Laplacian filter and re-normalization 
    13011348      DO jf = 1, nk 
     
    13141361      ENDIF 
    13151362#endif 
    1316     
     1363 
    13171364      ! AR(1) process and array swap 
    13181365      g2d = a*gb + b*g2d 
     
    13601407        ENDDO 
    13611408      ENDIF 
    1362        
     1409 
    13631410      ! Bound 
    13641411      IF( nn_sppt_step_bound .EQ. 2 ) THEN 
     
    14821529 
    14831530#ifdef NEMO_V34 
    1484       REWIND( numnam )             
     1531      REWIND( numnam ) 
    14851532      READ  ( numnam, namstopack ) 
    14861533#else 
    1487       REWIND( numnam_ref )  
     1534      REWIND( numnam_ref ) 
    14881535      READ  ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901) 
    14891536901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp ) 
    14901537 
    1491       REWIND( numnam_cfg )   
     1538      REWIND( numnam_cfg ) 
    14921539      READ  ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 ) 
    14931540902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp ) 
     
    15681615         WRITE(numout,*) 
    15691616         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       
     1617         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_spp_stdev :', rn_spp_stdev 
    15711618         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_spp_tau   :', rn_spp_tau 
    15721619         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    
     1620         WRITE(numout,*) '       SPP for bottom friction coeff                       nn_spp_bfr   :', nn_spp_bfr 
     1621         WRITE(numout,*) '                                            STDEV          rn_bfr_sd    :', rn_bfr_sd 
     1622         WRITE(numout,*) '       SPP for SST relaxation  coeff                       nn_spp_dqdt  :', nn_spp_dqdt 
     1623         WRITE(numout,*) '                                            STDEV          rn_dqdt_sd   :', rn_dqdt_sd 
     1624         WRITE(numout,*) '       SPP for SSS relaxation  coeff                       nn_spp_dedt  :', nn_spp_dedt 
     1625         WRITE(numout,*) '                                            STDEV          rn_dedt_sd   :', rn_dedt_sd 
     1626         WRITE(numout,*) '       SPP for vertical tra mixing coeff (only TKE, GLS)   nn_spp_avt   :', nn_spp_avt 
     1627         WRITE(numout,*) '                                            STDEV          rn_avt_sd    :', rn_avt_sd 
     1628         WRITE(numout,*) '       SPP for vertical dyn mixing coeff (only TKE, GLS)   nn_spp_avm   :', nn_spp_avm 
     1629         WRITE(numout,*) '                                            STDEV          rn_avm_sd    :', rn_avm_sd 
    15831630         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    
     1631         WRITE(numout,*) '                                            STDEV          rn_qsi0_sd   :', rn_qsi0_sd 
    15851632         WRITE(numout,*) '       SPP for horiz. diffusivity  U                       nn_spp_ahtu  :', nn_spp_ahtu 
    1586          WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd    
     1633         WRITE(numout,*) '                                            STDEV          rn_ahtu_sd   :', rn_ahtu_sd 
    15871634         WRITE(numout,*) '       SPP for horiz. diffusivity  V                       nn_spp_ahtv  :', nn_spp_ahtv 
    1588          WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd    
     1635         WRITE(numout,*) '                                            STDEV          rn_ahtv_sd   :', rn_ahtv_sd 
    15891636         WRITE(numout,*) '       SPP for horiz. diffusivity  W                       nn_spp_ahtw  :', nn_spp_ahtw 
    1590          WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd    
     1637         WRITE(numout,*) '                                            STDEV          rn_ahtw_sd   :', rn_ahtw_sd 
    15911638         WRITE(numout,*) '       SPP for horiz. diffusivity  T                       nn_spp_ahtt  :', nn_spp_ahtt 
    1592          WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd    
     1639         WRITE(numout,*) '                                            STDEV          rn_ahtt_sd   :', rn_ahtt_sd 
    15931640         WRITE(numout,*) '       SPP for horiz. viscosity (1/3)                      nn_spp_ahm1  :', nn_spp_ahm1 
    1594          WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd    
     1641         WRITE(numout,*) '                                            STDEV          rn_ahm1_sd   :', rn_ahm1_sd 
    15951642         WRITE(numout,*) '       SPP for horiz. viscosity (2/4)                      nn_spp_ahm2  :', nn_spp_ahm2 
    1596          WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd    
     1643         WRITE(numout,*) '                                            STDEV          rn_ahm2_sd   :', rn_ahm2_sd 
    15971644         WRITE(numout,*) '       SPP for relative wind factor                        nn_spp_relw  :', nn_spp_relw 
    15981645         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    
     1646         WRITE(numout,*) '                                            STDEV          rn_relw_sd   :', rn_relw_sd 
    16001647         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    
     1648         WRITE(numout,*) '                                            STDEV          rn_arnf_sd   :', rn_arnf_sd 
    16021649         WRITE(numout,*) '       SPP for geothermal heating                          nn_spp_geot  :', nn_spp_geot 
    1603          WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd    
     1650         WRITE(numout,*) '                                            STDEV          rn_geot_sd   :', rn_geot_sd 
    16041651         WRITE(numout,*) '       SPP for enhanced vertical diffusion                 nn_spp_aevd  :', nn_spp_aevd 
    1605          WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd    
     1652         WRITE(numout,*) '                                            STDEV          rn_aevd_sd   :', rn_aevd_sd 
    16061653         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    
     1654         WRITE(numout,*) '                                            STDEV          rn_tkelc_sd  :', rn_tkelc_sd 
    16081655         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    
     1656         WRITE(numout,*) '                                            STDEV          rn_tkedf_sd  :', rn_tkedf_sd 
    16101657         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    
     1658         WRITE(numout,*) '                                            STDEV          rn_tkeds_sd  :', rn_tkeds_sd 
    16121659         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    
     1660         WRITE(numout,*) '                                            STDEV          rn_tkebb_sd  :', rn_tkebb_sd 
    16141661         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    
     1662         WRITE(numout,*) '                                            STDEV          rn_tkefr_sd  :', rn_tkefr_sd 
    16161663         WRITE(numout,*) '       SPP for BBL U  diffusivity                          nn_spp_ahubbl:', nn_spp_ahubbl 
    16171664         WRITE(numout,*) '                                            STDEV          rn_ahubbl_sd :', rn_ahubbl_sd 
     
    16261673         WRITE(numout,*) 
    16271674         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     
     1675         WRITE(numout,*) '       SKEB switch                                         ln_skeb      :', ln_skeb 
     1676         WRITE(numout,*) '       SKEB ratio of backscattered energy                  rn_skeb      :', rn_skeb 
    16301677         WRITE(numout,*) '       Frequency update for dissipation mask               nn_dcom_freq :', nn_dcom_freq 
    16311678         WRITE(numout,*) '       Numerical dissipation factor (resolut. dependent)   rn_kh        :', rn_kh 
    16321679         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       
     1680         WRITE(numout,*) '       Standard deviation of random generator (AR1 field)  rn_skeb_stdev:', rn_skeb_stdev 
    16341681         WRITE(numout,*) '       Decorr. time scale                     (AR1 field)  rn_skeb_tau  :', rn_skeb_tau 
    16351682         WRITE(numout,*) '       Option of convection energy dissipation             nn_dconv     :', nn_dconv 
     
    17521799 
    17531800      ! Find filter attenuation factor 
    1754     
     1801 
    17551802      flt_fac = sppt_fltfac( sppt_filter_pass ) 
    17561803      rdt_sppt = nn_rndm_freq * rn_rdt 
    1757     
     1804 
    17581805      IF( ln_sppt_taumap ) THEN 
    17591806         CALL iom_open ( 'sppt_tau_map', inum ) 
     
    17981845      gauss_b = 0._wp 
    17991846      ! Weigths 
    1800       gauss_w(:)    = 1.0_wp  
     1847      gauss_w(:)    = 1.0_wp 
    18011848      IF( nn_vwei .eq. 1 ) THEN 
    18021849        gauss_w(1)    = 0.0_wp 
     
    18611908      IF(lwp .and. ln_stopack_diags) & 
    18621909      CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    1863     
     1910 
    18641911   END SUBROUTINE stopack_init 
    18651912   ! 
     
    18741921      INTEGER :: id1, jseed 
    18751922      CHARACTER(LEN=10)   ::   clseed='spsd0_0000' 
    1876       INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type  
    1877       INTEGER(KIND=8)     ::   ivals(8)  
     1923      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type 
     1924      INTEGER(KIND=8)     ::   ivals(8) 
    18781925      REAL(wp)            ::   zrseed4(4)           ! RNG seeds in integer type 
    18791926      REAL(wp)            ::   zrseed2d(jpi,jpj) 
     
    19832030      !!--------------------------------------------------------------------- 
    19842031      ! 
    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) , &  
     2032      ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,& 
     2033      gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , & 
    19872034      spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,& 
    19882035      gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),& 
     
    22082255      IF ( ln_dpsiv ) THEN 
    22092256       DO jp=1,jpni-1 
    2210          IF( jpri == jp ) THEN ! SEND TO EAST  
     2257         IF( jpri == jp ) THEN ! SEND TO EAST 
    22112258          zwrk(1:jpj) = dpsiv(jpi-1,:) 
    22122259          tag=2000+narea 
     
    22682315   REAL :: ds,dt,dtot,kh2 
    22692316   INTEGER :: ji,jj,jk 
    2270     
     2317 
    22712318   IF ( mt .eq. nit000 ) THEN 
    22722319        ALLOCATE ( dnum(jpi,jpj,jpk) ) 
     
    22872334          dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))*vmask(ji,jj,jk)*vmask(ji-1,jj,jk)/ e2v(ji,jj) + & 
    22882335               (un(ji,jj,jk)-un(ji,jj-1,jk))*umask(ji,jj,jk)*umask(ji,jj-1,jk)/ e1u(ji,jj) 
    2289    
     2336 
    22902337          dtot = sqrt( ds*ds + dt*dt ) * tmask(ji,jj,jk) 
    22912338          dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(ji,jj)*e2t(ji,jj) 
     
    22932340      ENDDO 
    22942341     ENDDO 
    2295     
     2342 
    22962343     CALL lbc_lnk(dnum,'T',1._wp) 
    22972344 
    22982345   ENDIF 
    22992346 
    2300    END SUBROUTINE  
     2347   END SUBROUTINE 
    23012348 
    23022349   SUBROUTINE skeb_dcon ( mt ) 
     
    23292376 
    23302377           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) )  
     2378              & / ( rau0 * fse3w(ji,jj,jk) ) 
    23322379 
    23332380           dcon(ji,jj,jk) = kc2*zz*e1t(ji,jj)*e2t(ji,jj)*rau0 / fse3w(ji,jj,jk) 
     
    23782425     IF(ln_skeb_own_gauss) THEN 
    23792426       DO jk=1,jpkm1 
    2380          psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2427         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:) 
    23812428       ENDDO 
    23822429     ELSE 
     
    24072454     IF(ln_skeb_own_gauss) THEN 
    24082455       DO jk=1,jpkm1 
    2409          psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2456         psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
    24102457       ENDDO 
    24112458     ELSE 
     
    24402487   IF(ln_skeb_own_gauss) THEN 
    24412488     DO jk=1,jpkm1 
    2442        psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)  
     2489       psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:) 
    24432490     ENDDO 
    24442491   ELSE 
  • branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r11442 r13088  
    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_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/step.F90

    r11442 r13088  
    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.