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 15603 – NEMO

Changeset 15603


Ignore:
Timestamp:
2021-12-16T10:39:55+01:00 (2 years ago)
Author:
mattmartin
Message:

Updated NEMO branch for coupled NWP at GO6 to include stochastic model perturbations.
For more info see ticket: ​https://code.metoffice.gov.uk/trac/nwpscience/ticket/1125.

Location:
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM
Files:
19 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/CONFIG/SHARED/field_def.xml

    r9163 r15603  
    125125         <field id="Age"        long_name="Sea Water Age Since Surface contact" standard_name="sea_water_age_since_surface_contact"  unit="yr" grid_ref="grid_T_3D"   /> 
    126126         <field id="Age_e3t"    long_name="Age * e3t"                                 unit="yr * m"  > Age * e3t </field > 
     127 
     128         <!-- variables available with ln_stopack set to .true. -->          
     129         <field id="sppt_ran" long_name="Random field for SPPT" standard_name="sppt_random_field" unit="0-1"     /> 
     130         <field id="skeb_ran" long_name="Random field for SKEB" standard_name="skeb_random_field" unit="0-1"     /> 
     131         <field id="spp_ran" long_name="Random field for SPP" standard_name="spp_random_field" unit="0-1"     /> 
     132         <field id="sppt_ar1" long_name="Perturbation field for SPPT" standard_name="sppt_perturb_field" unit="0-1" grid_ref="grid_T_3D"/> 
     133         <field id="spp_ar1" long_name="Perturbation field for SPP" standard_name="spp_perturb_field" unit="0-1" /> 
     134         <field id="skeb_ar1" long_name="Perturbation field for SKEB" standard_name="skeb_perturb_field" unit="0-1" /> 
     135 
     136         <!-- The following ones only when ln_stopack_diags is TRUE --> 
     137 
     138         <field id="spp_par01" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     139         <field id="spp_par02" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     140         <field id="spp_par03" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     141         <field id="spp_par04" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     142         <field id="spp_par05" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     143         <field id="spp_par06" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     144         <field id="spp_par07" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     145         <field id="spp_par08" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     146         <field id="spp_par09" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     147         <field id="spp_par10" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     148         <field id="spp_par11" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     149         <field id="spp_par12" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     150         <field id="spp_par13" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     151         <field id="spp_par14" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     152         <field id="spp_par15" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     153         <field id="spp_par16" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     154         <field id="spp_par17" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     155         <field id="spp_par18" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     156         <field id="spp_par19" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     157         <field id="spp_par20" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     158         <field id="spp_par21" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     159         <field id="spp_par22" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     160         <field id="spp_par23" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     161         <field id="spp_par24" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     162         <field id="spp_par25" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     163         <field id="spp_par26" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     164         <field id="spp_par27" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     165         <field id="spp_par28" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     166         <field id="spp_par29" long_name="Sample parameter SPP" standard_name="spp_pert_paramter" unit="0-1"     /> 
     167 
     168         <field id="ustar_skeb" long_name="Sea Water X Velocity Perturbation" standard_name="sea_water_x_velocity" unit="m/s" grid_ref="grid_U_3D" /> 
     169         <field id="vstar_skeb" long_name="Sea Water Y Velocity Perturbation" standard_name="sea_water_y_velocity" unit="m/s" grid_ref="grid_V_3D" /> 
     170 
    127171      </field_group> 
    128172 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r12555 r15603  
    3030   USE wrk_nemo        ! Memory Allocation 
    3131   USE timing          ! Timing 
     32   USE stopack 
    3233 
    3334   IMPLICIT NONE 
    … …  
    3839 
    3940   INTEGER ::   nldf = -2   ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals) 
     41 
     42#if defined key_dynldf_c3d 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahm10, ahm20, ahm30, ahm40 
     44#endif 
     45#if defined key_dynldf_c2d 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: ahm10, ahm20, ahm30, ahm40 
     47#endif 
    4048 
    4149   !! * Substitutions 
    … …  
    6775         ztrdv(:,:,:) = va(:,:,:)  
    6876      ENDIF 
     77 
     78#if defined key_dynldf_c3d 
     79         IF( kt == nit000 .AND. ln_stopack .AND. & 
     80            & ( (nn_spp_ahm1 + nn_spp_ahm2) > 0 ) ) THEN 
     81             ALLOCATE ( ahm10(jpi,jpj,jpk), ahm20(jpi,jpj,jpk) ) 
     82             ALLOCATE ( ahm30(jpi,jpj,jpk), ahm40(jpi,jpj,jpk) ) 
     83             ahm10 = ahm1 
     84             ahm20 = ahm2 
     85             ahm30 = ahm3 
     86             ahm40 = ahm4 
     87         ENDIF 
     88#endif 
     89#if defined key_dynldf_c2d 
     90         IF( kt == nit000 .AND. ln_stopack .AND. & 
     91            & ( (nn_spp_ahm1 + nn_spp_ahm2) > 0 ) ) THEN 
     92             ALLOCATE ( ahm10(jpi,jpj), ahm20(jpi,jpj) ) 
     93             ALLOCATE ( ahm30(jpi,jpj), ahm40(jpi,jpj) ) 
     94             ahm10 = ahm1 
     95             ahm20 = ahm2 
     96             ahm30 = ahm3 
     97             ahm40 = ahm4 
     98         ENDIF 
     99#endif 
     100 
     101#if defined key_traldf_c3d || defined key_traldf_c2d 
     102         IF( ln_stopack ) THEN 
     103            IF( nn_spp_ahm1 > 0 ) THEN 
     104               IF( ln_dynldf_lap ) THEN 
     105                  ahm1 = ahm10 
     106                  CALL spp_ahm(kt, ahm1, nn_spp_ahm1, rn_ahm1_sd, jk_spp_ahm1) 
     107               ENDIF 
     108               IF( ln_dynldf_bilap ) THEN 
     109                  ahm3 = ahm30 
     110                  CALL spp_ahm(kt, ahm3, nn_spp_ahm1, rn_ahm1_sd, jk_spp_ahm3) 
     111               ENDIF 
     112            ENDIF 
     113            IF( nn_spp_ahm2 > 0 ) THEN 
     114               IF( ln_dynldf_lap ) THEN 
     115               ahm2 = ahm20 
     116               CALL spp_ahm(kt, ahm2, nn_spp_ahm2, rn_ahm2_sd, jk_spp_ahm2) 
     117               ENDIF 
     118               IF( ln_dynldf_bilap ) THEN 
     119                  ahm4 = ahm40 
     120                  CALL spp_ahm(kt, ahm4, nn_spp_ahm2, rn_ahm2_sd, jk_spp_ahm4) 
     121               ENDIF 
     122            ENDIF 
     123         ENDIF 
     124#endif 
    69125 
    70126      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r12555 r15603  
    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) 
     24   USE stopack 
    2425 
    2526   IMPLICIT NONE 
    … …  
    3031 
    3132   INTEGER  ::   albd_init = 0      !: control flag for initialization 
    32    
     33 
    3334   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude 
    3435   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
    … …  
    3637   REAL(wp) ::   c2       = 0.10    !  "        " 
    3738   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
    38   
     39 
    3940   !                             !!* namelist namsbc_alb 
    4041   INTEGER  ::   nn_ice_alb 
    … …  
    5152      !!---------------------------------------------------------------------- 
    5253      !!               ***  ROUTINE albedo_ice  *** 
    53       !!           
    54       !! ** Purpose :   Computation of the albedo of the snow/ice system  
    55       !!        
     54      !! 
     55      !! ** Purpose :   Computation of the albedo of the snow/ice system 
     56      !! 
    5657      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
    5758      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
    … …  
    7273      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
    7374      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
    74       !!                  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 
    7576      !!                     under melting conditions than under freezing conditions 
    76       !!                  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 
    7778      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 
    7879      !! 
    7980      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 
    8081      !!                Brandt et al. 2005, J. Climate, vol 18 
    81       !!                Grenfell & Perovich 2004, JGR, vol 109  
     82      !!                Grenfell & Perovich 2004, JGR, vol 109 
    8283      !!---------------------------------------------------------------------- 
    8384      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin) 
    … …  
    9697 
    9798      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    98        
     99 
    99100      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    100101 
    101       IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    102  
    103        
     102      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
     103 
     104 
    104105      SELECT CASE ( nn_ice_alb ) 
    105106 
    … …  
    108109      !------------------------------------------ 
    109110      CASE( 0 ) 
    110         
     111 
    111112         ralb_sf = 0.80       ! dry snow 
    112113         ralb_sm = 0.65       ! melting snow 
    113114         ralb_if = 0.72       ! bare frozen ice 
    114          ralb_im = rn_albice  ! bare puddled ice  
    115           
     115         ralb_im = rn_albice  ! bare puddled ice 
     116 
    116117         !  Computation of ice albedo (free of snow) 
    117118         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
    118119         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if 
    119120         END  WHERE 
    120        
     121 
    121122         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    122123         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
    … …  
    126127         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
    127128         END WHERE 
    128       
     129 
    129130         DO jl = 1, ijpl 
    130131            DO jj = 1, jpj 
    … …  
    132133                  ! freezing snow 
    133134                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
    134                   !                                        !  freezing snow         
     135                  !                                        !  freezing snow 
    135136                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    136137                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
    137138                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   & 
    138                      &        +         zswitch   * ralb_sf   
     139                     &        +         zswitch   * ralb_sf 
    139140 
    140141                  ! melting snow 
    … …  
    142143                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
    143144                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   & 
    144                       &     +         zswitch   *   ralb_sm  
     145                      &     +         zswitch   *   ralb_sm 
    145146                  ! 
    146147                  ! snow albedo 
    147                   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 ) ) 
    148149                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    149                 
     150 
    150151                  ! Ice/snow albedo 
    151152                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    … …  
    154155               END DO 
    155156            END DO 
     157 
     158#if defined key_traldf_c2d || key_traldf_c3d 
     159            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
     160               & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 
     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 
    156166         END DO 
    157167 
    … …  
    161171      !  New parameterization (2016) 
    162172      !------------------------------------------ 
    163       CASE( 1 )  
     173      CASE( 1 ) 
    164174 
    165175         ralb_im = rn_albice  ! bare puddled ice 
    … …  
    176186!         ralb_sm = 0.82      ! melting snow 
    177187!         ralb_if = 0.54      ! bare frozen ice 
    178 !  
     188! 
    179189         !  Computation of ice albedo (free of snow) 
    180          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     190         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
    181191         z1_c2 = 1. / 0.05 
    182192         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im 
    183193         ELSE WHERE                                              ;   zalb = ralb_if 
    184194         END  WHERE 
    185           
     195 
    186196         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    187197         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
    … …  
    200210 
    201211                   ! snow albedo 
    202                   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 ) ) 
    203213                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    204214 
    205                   ! Ice/snow albedo    
     215                  ! Ice/snow albedo 
    206216                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    207217                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
    … …  
    209219              END DO 
    210220            END DO 
     221 
     222#if defined key_traldf_c2d || key_traldf_c3d 
     223            IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 
     224               & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 
     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 
    211230         END DO 
    212231         ! Effect of the clouds (2d order polynomial) 
    213          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 ); 
    214233 
    215234      END SELECT 
    216        
     235 
    217236      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    218237      ! 
    … …  
    223242      !!---------------------------------------------------------------------- 
    224243      !!               ***  ROUTINE albedo_oce  *** 
    225       !!  
     244      !! 
    226245      !! ** Purpose :   Computation of the albedo of the ocean 
    227246      !!---------------------------------------------------------------------- 
    … …  
    229248      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky 
    230249      !! 
    231       REAL(wp) :: zcoef  
     250      REAL(wp) :: zcoef 
    232251      !!---------------------------------------------------------------------- 
    233252      ! 
    234253      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982 
    235       pa_oce_cs(:,:) = zcoef  
     254      pa_oce_cs(:,:) = zcoef 
    236255      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972 
    237256      ! 
    … …  
    248267      !!---------------------------------------------------------------------- 
    249268      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    250       NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
     269      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
    251270      !!---------------------------------------------------------------------- 
    252271      ! 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r12555 r15603  
    5151   USE par_ice_2 
    5252#endif 
     53   USE stopack 
    5354 
    5455   IMPLICIT NONE 
    … …  
    8990   REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem) 
    9091   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
     92   REAL(wp), ALLOCATABLE, SAVE ::   rn_vfac0(:,:) ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
    9193   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9294   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     95   REAL(wp), PUBLIC :: rn_sfac ! multiplication factor for snow precipitation over sea-ice 
    9396 
    9497   !! * Substitutions 
    … …  
    151154         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    152155         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    153          &                  sn_tdif, rn_zqt,  rn_zu 
     156         &                  sn_tdif, rn_zqt,  rn_zu, rn_sfac 
    154157      !!--------------------------------------------------------------------- 
    155158      ! 
    … …  
    158161         !                                      ! ====================== ! 
    159162         ! 
     163         rn_sfac = 1._wp       ! Default to one if missing from namelist 
    160164         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
    161165         READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) 
    … …  
    196200         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    197201         ! 
    198       ENDIF 
     202         ALLOCATE ( rn_vfac0(jpi,jpj) ) 
     203         rn_vfac0(:,:) = rn_vfac 
     204         ! 
     205      ENDIF 
     206 
     207#if defined key_traldf_c2d || key_traldf_c3d 
     208      IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 
     209         rn_vfac0(:,:) = rn_vfac 
     210         CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 
     211      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 
    199217 
    200218      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
    … …  
    205223#if defined key_cice 
    206224      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    207          qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)  
     225         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
    208226         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    209227         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
    210228         ENDIF 
    211          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)          
     229         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    212230         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
    213231         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    … …  
    219237      ! 
    220238   END SUBROUTINE sbc_blk_core 
    221     
    222     
     239 
     240 
    223241   SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 
    224242      !!--------------------------------------------------------------------- 
    … …  
    230248      !! ** Method  :   CORE bulk formulea for the ocean using atmospheric 
    231249      !!      fields read in sbc_read 
    232       !!  
     250      !! 
    233251      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
    234252      !!              - vtau    : j-component of the stress at V-point  (N/m2) 
    … …  
    268286      ! local scalars ( place there for vector optimisation purposes) 
    269287      zcoef_qsatw = 0.98 * 640380. / rhoa 
    270        
     288 
    271289      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    272290 
    … …  
    276294 
    277295      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    278       zwnd_i(:,:) = 0.e0   
     296      zwnd_i(:,:) = 0.e0 
    279297      zwnd_j(:,:) = 0.e0 
    280298#if defined key_cyclone 
    … …  
    289307      DO jj = 2, jpjm1 
    290308         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    291             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    292             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     309            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     310            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    293311         END DO 
    294312      END DO 
    … …  
    320338      CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
    321339         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    322      
     340 
    323341      ! ... tau module, i and j component 
    324342      DO jj = 1, jpj 
    … …  
    351369      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    352370 
    353      
     371 
    354372      !  Turbulent fluxes over ocean 
    355373      ! ----------------------------- 
    … …  
    376394         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce_core: zst    : ') 
    377395      ENDIF 
    378         
     396 
    379397      ! ----------------------------------------------------------------------------- ! 
    380398      !     III    Total FLUXES                                                       ! 
    … …  
    384402         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    385403      ! 
    386       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
     404      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    387405         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    388406         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    … …  
    425443      ! 
    426444   END SUBROUTINE blk_oce_core 
    427   
    428     
     445 
     446 
    429447#if defined key_lim2 || defined key_lim3 
    430448   SUBROUTINE blk_ice_core_tau 
    … …  
    465483            DO ji = 2, jpim1   ! B grid : NO vector opt 
    466484               ! ... scalar wind at I-point (fld being at T-point) 
    467                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    468                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    469                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    470                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     485               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)    & 
     486                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) & 
     487                  &      - rn_vfac0(ji,jj) * u_ice(ji,jj) 
     488               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)    & 
     489                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) & 
     490                  &      - rn_vfac0(ji,jj) * v_ice(ji,jj) 
    471491               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    472492               ! ... ice stress at I-point 
    … …  
    474494               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    475495               ! ... scalar wind at T-point (fld being at T-point) 
    476                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
    477                   &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
    478                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    479                   &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     496               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1)                                         & 
     497                  &      - rn_vfac0(ji,jj) * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     498                  &                                  + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     499               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1)                                         & 
     500                  &      - rn_vfac0(ji,jj) * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     501                  &                                  + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    480502               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    481503            END DO 
    … …  
    488510         DO jj = 2, jpj 
    489511            DO ji = fs_2, jpi   ! vect. opt. 
    490                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    491                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     512               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     513               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac0(ji,jj) * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    492514               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    493515            END DO 
    … …  
    495517         DO jj = 2, jpjm1 
    496518            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    497                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    498                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    499                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    500                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
     519               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )             & 
     520                  &              * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) & 
     521                  &                  - rn_vfac0(ji,jj) * u_ice(ji,jj) ) 
     522               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )             & 
     523                  &              * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) & 
     524                  &                  - rn_vfac0(ji,jj) * v_ice(ji,jj) ) 
    501525            END DO 
    502526         END DO 
    … …  
    513537 
    514538      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
    515        
     539 
    516540   END SUBROUTINE blk_ice_core_tau 
    517541 
    … …  
    526550      !!                between atmosphere and sea-ice using CORE bulk 
    527551      !!                formulea, ice variables and read atmmospheric fields. 
    528       !!  
     552      !! 
    529553      !! caution : the net upward water flux has with mm/day unit 
    530554      !!--------------------------------------------------------------------- 
    … …  
    546570      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
    547571      ! 
    548       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 ) 
    549573 
    550574      ! local scalars ( place there for vector optimisation purposes) 
    … …  
    569593               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    570594               ! lw sensitivity 
    571                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     595               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
    572596 
    573597               ! ----------------------------! 
    … …  
    579603               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    580604               ! Latent Heat 
    581                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)   & 
    582606                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    583607              ! Latent heat sensitivity for ice (Dqla/Dt) 
    … …  
    610634 
    611635#if defined  key_lim3 
    612       CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     636      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
    613637 
    614638      ! --- evaporation --- ! 
    … …  
    620644      ! --- evaporation minus precipitation --- ! 
    621645      zsnw(:,:) = 0._wp 
    622       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 
    623647      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    624648      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
    … …  
    643667      DO jl = 1, jpl 
    644668         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
    645                                    ! But we do not have Tice => consider it at 0°C => evap=0  
     669                                   ! But we do not have Tice => consider it at 0 degC => evap=0 
    646670      END DO 
    647671 
    648       CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     672      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
    649673#endif 
    650674 
    … …  
    670694      ! 
    671695      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
    672        
     696 
    673697   END SUBROUTINE blk_ice_core_flx 
    674698#endif 
    … …  
    683707      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    684708      !! 
    685       !! ** Method : Monin Obukhov Similarity Theory  
     709      !! ** Method : Monin Obukhov Similarity Theory 
    686710      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    687711      !! 
    … …  
    693717      !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    694718      !!       rather than temperature difference only... 
    695       !!    - added function "cd_neutral_10m" that uses the improved parametrization of  
     719      !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    696720      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    697721      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
    … …  
    728752 
    729753      IF( nn_timing == 1 )  CALL timing_start('turb_core_2z') 
    730      
     754 
    731755      CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
    732756      CALL wrk_alloc( jpi,jpj, zeta_u, stab ) 
    … …  
    740764      U_zu = MAX( 0.5 , dU )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    741765 
    742       !! First guess of stability:  
     766      !! First guess of stability: 
    743767      ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 
    744768      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    … …  
    754778      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    755779      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    756      
     780 
    757781      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    758782      Cd = ztmp0   ;   Ce = Ce_n10   ;   Ch = Ch_n10   ;   sqrt_Cd = sqrt_Cd_n10 
    … …  
    765789         ! 
    766790         ztmp1 = T_zu - sst   ! Updating air/sea differences 
    767          ztmp2 = q_zu - q_sat  
     791         ztmp2 = q_zu - q_sat 
    768792 
    769793         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    770794         ztmp1  = Ch/sqrt_Cd*ztmp1    ! theta* 
    771795         ztmp2  = Ce/sqrt_Cd*ztmp2    ! q* 
    772         
     796 
    773797         ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 
    774798 
    775799         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    776          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) 
    777801         !                                                                     ( Cd*U_zu*U_zu is U*^2 at zu) 
    778802 
    … …  
    781805         zpsi_h_u = psi_h( zeta_u ) 
    782806         zpsi_m_u = psi_m( zeta_u ) 
    783         
     807 
    784808         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 
    785809         IF ( .NOT. l_zt_equal_zu ) THEN 
    … …  
    790814            q_zu = max(0., q_zu) 
    791815         END IF 
    792         
     816 
    793817         IF( ln_cdgw ) THEN      ! surface wave case 
    794             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u )  
     818            sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
    795819            Cd      = sqrt_Cd * sqrt_Cd 
    796820         ELSE 
    … …  
    802826           ztmp0 = cd_neutral_10m(ztmp0)                                                 ! Cd_n10 
    803827           sqrt_Cd_n10 = sqrt(ztmp0) 
    804         
     828 
    805829           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                     ! L&Y 2004 eq. (6b) 
    806830           stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    … …  
    809833           !! Update of transfer coefficients: 
    810834           ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u)   ! L&Y 2004 eq. (10a) 
    811            Cd      = ztmp0 / ( ztmp1*ztmp1 )    
     835           Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    812836           sqrt_Cd = SQRT( Cd ) 
    813837         ENDIF 
    … …  
    815839         ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    816840         ztmp2 = sqrt_Cd / sqrt_Cd_n10 
    817          ztmp1 = 1. + Ch_n10*ztmp0                
     841         ztmp1 = 1. + Ch_n10*ztmp0 
    818842         Ch  = Ch_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    819843         ! 
    820          ztmp1 = 1. + Ce_n10*ztmp0                
     844         ztmp1 = 1. + Ce_n10*ztmp0 
    821845         Ce  = Ce_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    822846         ! 
    … …  
    836860   FUNCTION cd_neutral_10m( zw10 ) 
    837861      !!---------------------------------------------------------------------- 
    838       !! Estimate of the neutral drag coefficient at 10m as a function  
     862      !! Estimate of the neutral drag coefficient at 10m as a function 
    839863      !! of neutral wind  speed at 10m 
    840864      !! 
    … …  
    842866      !! 
    843867      !! Author: L. Brodeau, june 2014 
    844       !!----------------------------------------------------------------------     
     868      !!---------------------------------------------------------------------- 
    845869      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   zw10           ! scalar wind speed at 10m (m/s) 
    846870      REAL(wp), DIMENSION(jpi,jpj)             ::   cd_neutral_10m 
    847871      ! 
    848872      REAL(wp), DIMENSION(:,:), POINTER ::   rgt33 
    849       !!----------------------------------------------------------------------     
     873      !!---------------------------------------------------------------------- 
    850874      ! 
    851875      CALL wrk_alloc( jpi,jpj, rgt33 ) 
    852876      ! 
    853877      !! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    854       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 
    855879      cd_neutral_10m = 1.e-3 * ( & 
    856880         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r12555 r15603  
    5454   REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s] 
    5555   REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff 
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_avt_rnf0 
    5657 
    5758   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis 
    … …  
    159160            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 
    160161         ENDIF 
    161   
     162 
    162163         IF(lwp .AND. lflush) CALL flush(numout) 
    163164 
    … …  
    465466         !                                      !    - mixed upstream-centered (ln_traadv_cen2=T) 
    466467         ! 
     468         ALLOCATE ( rn_avt_rnf0(jpi,jpj) ) 
     469         rn_avt_rnf0(:,:) = rn_avt_rnf 
     470         ! 
    467471         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
    468472            &                                              'be spread through depth by ln_rnf_depth'               ) 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r12555 r15603  
    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) 
     27   USE stopack 
     28   USE wrk_nemo       ! Memory Allocation 
    2729 
    2830   IMPLICIT NONE 
    … …  
    4042   REAL(wp)        ::   rn_dqdt         ! restoring factor on SST and SSS 
    4143   REAL(wp)        ::   rn_deds         ! restoring factor on SST and SSS 
    42    LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term  
     44   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
    4345   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day] 
    4446 
    … …  
    7577      REAL(wp) ::   zerp     ! local scalar for evaporation damping 
    7678      REAL(wp) ::   zqrp     ! local scalar for heat flux damping 
    77       REAL(wp) ::   zsrp     ! local scalar for unit conversion of rn_deds factor 
    7879      REAL(wp) ::   zerp_bnd ! local scalar for unit conversion of rn_epr_max factor 
     80      REAL(wp), POINTER, DIMENSION(:,:) :: rn_dqdt_s, zsrp 
    7981      INTEGER  ::   ierror   ! return error code 
    8082      !! 
    … …  
    9597            ! 
    9698            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term 
     99 
     100               CALL wrk_alloc( jpi, jpj, rn_dqdt_s ) 
     101               rn_dqdt_s(:,:) = rn_dqdt 
     102 
     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 
    97112               DO jj = 1, jpj 
    98113                  DO ji = 1, jpi 
    99                      zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
     114                     zqrp = rn_dqdt_s(ji,jj) * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
    100115                     qns(ji,jj) = qns(ji,jj) + zqrp 
    101116                     qrp(ji,jj) = zqrp 
    … …  
    103118               END DO 
    104119               CALL iom_put( "qrp", qrp )                             ! heat flux damping 
     120               CALL wrk_dealloc( jpi, jpj, rn_dqdt_s ) 
    105121            ENDIF 
    106122            ! 
    107123            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx)) 
    108                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
     124               CALL wrk_alloc( jpi, jpj, zsrp) 
     125               zsrp(:,:) = rn_deds 
     126#if defined key_traldf_c2d || key_traldf_c3d 
     127               IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 
     128                  & 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 
    109136!CDIR COLLAPSE 
    110137               DO jj = 1, jpj 
    111138                  DO ji = 1, jpi 
    112                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    113                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )  
     139                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     140                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 
    114141                     sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux 
    115142                     erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only) 
    … …  
    117144               END DO 
    118145               CALL iom_put( "erp", erp )                             ! freshwater flux damping 
     146               CALL wrk_dealloc( jpi,jpj, zsrp ) 
    119147               ! 
    120148            ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns) 
    121                zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
    122                zerp_bnd = rn_sssr_bnd / rday                          !       -              -     
     149               CALL wrk_alloc( jpi, jpj, zsrp) 
     150               zsrp(:,:) = rn_deds 
     151#if defined key_traldf_c2d || key_traldf_c3d 
     152               IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 
     153                  & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 
     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                          !       -              - 
    123160!CDIR COLLAPSE 
    124161               DO jj = 1, jpj 
    125                   DO ji = 1, jpi                             
    126                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
     162                  DO ji = 1, jpi 
     163                     zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    127164                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    128165                        &        / MAX(  sss_m(ji,jj), 1.e-20   ) 
    … …  
    134171               END DO 
    135172               CALL iom_put( "erp", erp )                             ! freshwater flux damping 
     173               CALL wrk_dealloc( jpi,jpj,zsrp ) 
    136174            ENDIF 
    137175            ! 
    … …  
    144182   END SUBROUTINE sbc_ssr 
    145183 
    146   
     184 
    147185   SUBROUTINE sbc_ssr_init 
    148186      !!--------------------------------------------------------------------- 
    … …  
    167205      !!---------------------------------------------------------------------- 
    168206      ! 
    169   
    170       REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist :  
     207 
     208      REWIND( numnam_ref )              ! Namelist namsbc_ssr in reference namelist : 
    171209      READ  ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901) 
    172210901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp ) 
    … …  
    224262      ENDIF 
    225263      ! 
    226       !                            !* Initialize qrp and erp if no restoring  
     264      !                            !* Initialize qrp and erp if no restoring 
    227265      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp 
    228266      IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp 
    229267      ! 
    230268   END SUBROUTINE sbc_ssr_init 
    231        
     269 
    232270   !!====================================================================== 
    233271END MODULE sbcssr 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r12555 r15603  
    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 
    … …  
    2828   USE wrk_nemo        ! Memory Allocation 
    2929   USE timing          ! Timing 
     30   USE stopack 
    3031 
    3132   IMPLICIT NONE 
    … …  
    4142 
    4243   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend 
     44   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd1   ! geothermal heating trend 
    4345   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read) 
    44   
     46 
    4547   !! * Substitutions 
    4648#  include "domzgr_substitute.h90" 
    … …  
    5658      !!                  ***  ROUTINE tra_bbc  *** 
    5759      !! 
    58       !! ** Purpose :   Compute the bottom boundary contition on temperature  
    59       !!              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 
    6062      !!              general trend of temperature equations. 
    6163      !! 
    62       !! ** Method  :   The geothermal heat flux set to its constant value of  
     64      !! ** Method  :   The geothermal heat flux set to its constant value of 
    6365      !!              86.4 mW/m2 (Stein and Stein 1992, Huang 1999). 
    6466      !!       The temperature trend associated to this heat flux through the 
    … …  
    8991      ! 
    9092      !                             !  Add the geothermal heat flux trend on temperature 
     93 
     94#if defined key_traldf_c2d || key_traldf_c3d 
     95      IF( ln_stopack .AND. nn_spp_geot > 0) THEN 
     96          qgh_trd1(:,:) = qgh_trd0(:,:) 
     97          CALL spp_gen(kt, qgh_trd1, nn_spp_geot, rn_geot_sd, jk_spp_geot) 
     98      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 
    91106      DO jj = 2, jpjm1 
    92107         DO ji = 2, jpim1 
    93108            ik = mbkt(ji,jj) 
    94             zqgh_trd = qgh_trd0(ji,jj) / fse3t(ji,jj,ik) 
     109            zqgh_trd = qgh_trd1(ji,jj) / fse3t(ji,jj,ik) 
    95110            tsa(ji,jj,ik,jp_tem) = tsa(ji,jj,ik,jp_tem) + zqgh_trd 
    96111         END DO 
    … …  
    137152      CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files 
    138153      ! 
    139       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 
    140155      !!---------------------------------------------------------------------- 
    141156 
    … …  
    164179         ! 
    165180         ALLOCATE( qgh_trd0(jpi,jpj) )    ! allocation 
     181         ALLOCATE( qgh_trd1(jpi,jpj) )    ! allocation 
    166182         ! 
    167183         SELECT CASE ( nn_geoflx )        ! geothermal heat flux / (rauO * Cp) 
    … …  
    193209            ! 
    194210         END SELECT 
     211         qgh_trd1(:,:) = qgh_trd0(:,:) 
    195212         ! 
    196213      ELSE 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r12555 r15603  
    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) 
     41   USE stopack 
    4142 
    4243   IMPLICIT NONE 
    … …  
    6768   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM) 
    6869   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points 
     70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_1, ahv_bbl_1   ! diffusive bbl flux coefficients at u and v-points 
    6971   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM) 
    7072 
    … …  
    8688         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     & 
    8789         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          & 
     90         &      ahu_bbl_1(jpi,jpj) , ahv_bbl_1(jpi,jpj) ,                                          & 
    8891         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc ) 
    8992         ! 
    … …  
    195198      ALLOCATE(zptb(1:jpi, 1:jpj)) 
    196199      ! 
     200      ahu_bbl_1(:,:) = ahu_bbl(:,:) 
     201#if defined key_traldf_c2d || key_traldf_c3d 
     202      IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN 
     203          CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl ) 
     204      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 
     212      ahv_bbl_1(:,:) = ahv_bbl(:,:) 
     213#if defined key_traldf_c2d || key_traldf_c3d 
     214      IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN 
     215          CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl ) 
     216      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 
     223      ! 
    197224      DO jn = 1, kjpt                                     ! tracer loop 
    198225         !                                                ! =========== 
    … …  
    203230            END DO 
    204231         END DO 
    205          !                
     232         ! 
    206233         DO jj = 2, jpjm1                                    ! Compute the trend 
    207234            DO ji = 2, jpim1 
    … …  
    209236               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik) 
    210237               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
    211                   &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
    212                   &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
    213                   &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
    214                   &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
     238                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   & 
     239                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   & 
     240                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   & 
     241                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr 
    215242            END DO 
    216243         END DO 
    … …  
    415442                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
    416443                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
    417                   !                                                          ! 2*masked bottom density gradient  
     444                  !                                                          ! 2*masked bottom density gradient 
    418445                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
    419446                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
    … …  
    578605               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    579606            ENDIF 
    580             !      
     607            ! 
    581608            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN 
    582609               mgrhv(ji,jj) = INT(  SIGN( 1.e0, & 
    … …  
    598625      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1) 
    599626      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1) 
    600  
    601627 
    602628      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r12555 r15603  
    3232   USE wrk_nemo        ! Memory allocation 
    3333   USE timing          ! Timing 
     34   USE stopack 
    3435 
    3536   IMPLICIT NONE 
    … …  
    4344   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile 
    4445   !                                                               !  (key_traldf_ano only) 
     46#if defined key_traldf_c3d 
     47   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: ahtu0, ahtv0, ahtw0, ahtt0 
     48#endif 
     49#if defined key_traldf_c2d 
     50   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:  ) :: ahtu0, ahtv0, ahtw0, ahtt0 
     51#endif 
    4552 
    4653   !! * Substitutions 
    … …  
    7582         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    7683      ENDIF 
     84 
     85#if defined key_traldf_c3d 
     86         IF(  ( kt == nit000 ) .AND. & 
     87            & ( ln_stopack )   .AND. & 
     88            & ( ( nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt ) > 0 ) ) THEN 
     89             ALLOCATE ( ahtu0(jpi,jpj,jpk), ahtv0(jpi,jpj,jpk) ) 
     90             ALLOCATE ( ahtt0(jpi,jpj,jpk), ahtw0(jpi,jpj,jpk) ) 
     91             ahtu0 = ahtu 
     92             ahtv0 = ahtv 
     93             ahtw0 = ahtw 
     94             ahtt0 = ahtt 
     95         ENDIF 
     96#endif 
     97#if defined key_traldf_c2d 
     98         IF(  ( kt == nit000 ) .AND. & 
     99            & ( ln_stopack )   .AND. & 
     100            & ( ( nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt ) > 0 ) ) THEN 
     101             ALLOCATE ( ahtu0(jpi,jpj), ahtv0(jpi,jpj) ) 
     102             ALLOCATE ( ahtt0(jpi,jpj), ahtw0(jpi,jpj) ) 
     103             ahtu0 = ahtu 
     104             ahtv0 = ahtv 
     105             ahtw0 = ahtw 
     106             ahtt0 = ahtt 
     107         ENDIF 
     108#endif 
     109#if defined key_traldf_c3d || defined key_traldf_c2d 
     110         IF( ln_stopack .AND. ( nn_spp_ahtu > 0 ) ) THEN 
     111             ahtu = ahtu0 
     112             CALL spp_aht(kt, ahtu, nn_spp_ahtu, rn_ahtu_sd, jk_spp_ahtu) 
     113         ENDIF 
     114         IF( ln_stopack .AND. ( nn_spp_ahtv > 0 ) ) THEN 
     115             ahtv = ahtv0 
     116             CALL spp_aht(kt, ahtv, nn_spp_ahtv, rn_ahtv_sd, jk_spp_ahtv) 
     117         ENDIF 
     118         IF( ln_stopack .AND. ( nn_spp_ahtw > 0 ) ) THEN 
     119             ahtw = ahtw0 
     120             CALL spp_aht(kt, ahtw, nn_spp_ahtw, rn_ahtw_sd, jk_spp_ahtw) 
     121         ENDIF 
     122         IF( ln_stopack .AND. ( nn_spp_ahtt > 0 ) ) THEN 
     123             ahtt = ahtt0 
     124             CALL spp_aht(kt, ahtt, nn_spp_ahtt, rn_ahtt_sd, jk_spp_ahtt) 
     125         ENDIF 
     126#endif 
    77127 
    78128      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r12555 r15603  
    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   !!---------------------------------------------------------------------- 
    … …  
    3333   USE wrk_nemo       ! Memory Allocation 
    3434   USE timing         ! Timing 
     35   USE stopack 
    3536 
    3637   IMPLICIT NONE 
    … …  
    4243   !                                 !!* Namelist namtra_qsr: penetrative solar radiation 
    4344   LOGICAL , PUBLIC ::   ln_traqsr    !: light absorption (qsr) flag 
    44    LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag   
     45   LOGICAL , PUBLIC ::   ln_qsr_rgb   !: Red-Green-Blue light absorption flag 
    4546   LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag 
    4647   LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag 
    … …  
    5051   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    5152   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands) 
    52   
     53 
    5354   ! Module variables 
    54    REAL(wp) ::   xsi0r                           !: inverse of rn_si0 
     55   REAL(wp), ALLOCATABLE ::   xsi0r(:,:)         !: inverse of rn_si0 
    5556   REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
    5657   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    … …  
    7980      !!      Considering the 2 wavebands case: 
    8081      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
    81       !!         The temperature trend associated with the solar radiation penetration  
     82      !!         The temperature trend associated with the solar radiation penetration 
    8283      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 
    8384      !!         At the bottom, boudary condition for the radiation is no flux : 
    … …  
    8586      !!      in the last ocean level. 
    8687      !!         In z-coordinate case, the computation is only done down to the 
    87       !!      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 
    8889      !!      used for the computation are calculated one for once as they 
    8990      !!      depends on k only. 
    … …  
    105106      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         - 
    106107      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
    107       REAL(wp) ::   zlogc, zlogc2, zlogc3  
     108      REAL(wp) ::   zlogc, zlogc2, zlogc3 
    108109      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
    109110      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d 
    … …  
    112113      IF( nn_timing == 1 )  CALL timing_start('tra_qsr') 
    113114      ! 
    114       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    115       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 ) 
    116117      ! 
    117118      IF( kt == nit000 ) THEN 
    … …  
    124125 
    125126      IF( l_trdtra ) THEN      ! Save ta and sa trends 
    126          CALL wrk_alloc( jpi, jpj, jpk, ztrdt )  
     127         CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 
    127128         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    128129      ENDIF 
    … …  
    153154      !                                        Compute now qsr tracer content field 
    154155      !                                        ************************************ 
    155        
     156 
    156157      !                                           ! ============================================== ! 
    157158      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  ! 
    … …  
    162163         !                                        Add to the general trend 
    163164         DO jk = 1, jpkm1 
    164             DO jj = 2, jpjm1  
     165            DO jj = 2, jpjm1 
    165166               DO ji = fs_2, fs_jpim1   ! vector opt. 
    166167                  z1_e3t = zfact / fse3t(ji,jj,jk) 
    … …  
    183184         ENDIF 
    184185         !                                        ! ============================================== ! 
    185       ELSE                                        !  Ocean alone :  
     186      ELSE                                        !  Ocean alone : 
    186187         !                                        ! ============================================== ! 
    187188         ! 
    188          !                                                ! ------------------------- ! 
     189         ! 
     190#if defined key_traldf_c2d || key_traldf_c3d 
     191         IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN 
     192            xsi0r = rn_si0 
     193            CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 
     194            xsi0r = 1.e0 / xsi0r 
     195         ENDIF 
     196#else 
     197         IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) & 
     198            & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// & 
     199                             'key_traldf_c2d or key_traldf_c3d') 
     200#endif 
     201 
     202         !                                               ! ------------------------- ! 
    189203         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration ! 
    190204            !                                             ! ------------------------- ! 
    … …  
    196210                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step 
    197211                  DO jk = 1, nksr + 1 
    198                      zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)  
     212                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
    199213                  ENDDO 
    200214                  ! 
    … …  
    217231                        zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
    218232                        zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
    219                         zCze    = 1.12  * (zchl)**0.803  
     233                        zCze    = 1.12  * (zchl)**0.803 
    220234                        DO jk = 1, nksr + 1 
    221235                           zpsi = fsdept(ji,jj,jk) / zze 
    … …  
    228242               ELSE                              !* Variable ocean volume but constant chrlorophyll 
    229243                  DO jk = 1, nksr + 1 
    230                      zchl3d(:,:,jk) = 0.05  
     244                     zchl3d(:,:,jk) = 0.05 
    231245                  ENDDO 
    232246               ENDIF 
    … …  
    253267!CDIR NOVERRCHK 
    254268                  DO jj = 1, jpj 
    255 !CDIR NOVERRCHK    
    256                      DO ji = 1, jpi 
    257                         zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r    ) 
     269!CDIR NOVERRCHK 
     270                     DO ji = 1, jpi 
     271                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) 
    258272                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 
    259273                        zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 
    … …  
    286300                     END DO 
    287301                  END DO 
    288                   !  
     302                  ! 
    289303                  DO jj = 1, jpj 
    290304                     DO ji = 1, jpi 
    291                         zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r    ) 
     305                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r(ji,jj) ) 
    292306                        zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 
    293307                        zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
    294308                        zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
    295                         fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)  
     309                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2) 
    296310                     END DO 
    297311                  END DO 
    … …  
    314328            !                                             ! ------------------------- ! 
    315329            ! 
    316             IF( lk_vvl ) THEN                                  !* variable volume 
     330            IF( lk_vvl .OR. ( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) ) THEN        !* variable volume 
     331 
    317332               zz0   =        rn_abs   * r1_rau0_rcp 
    318333               zz1   = ( 1. - rn_abs ) * r1_rau0_rcp 
    319                DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m  
     334               DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m 
    320335                  DO jj = 1, jpj 
    321336                     DO ji = 1, jpi 
    322                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    323                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    324                         qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )  
     337                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
     338                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
     339                        qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 
    325340                     END DO 
    326341                  END DO 
    … …  
    330345                  DO jj = 1, jpj 
    331346                     DO ji = 1, jpi 
    332                         zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
    333                         zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
     347                        zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 
     348                        zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 
    334349                        fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 
    335350                     END DO 
    … …  
    340355                  DO jj = 2, jpjm1 
    341356                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    342                         ! (ISF) no light penetration below the ice shelves          
     357                        ! (ISF) no light penetration below the ice shelves 
    343358                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 
    344359                     END DO 
    … …  
    356371         !                                        Add to the general trend 
    357372         DO jk = 1, nksr 
    358             DO jj = 2, jpjm1  
     373            DO jj = 2, jpjm1 
    359374               DO ji = fs_2, fs_jpim1   ! vector opt. 
    360375                  z1_e3t = zfact / fse3t(ji,jj,jk) 
    … …  
    375390            IF(lflush) CALL flush(numout) 
    376391         ENDIF 
    377          IF(nn_timing == 2)  CALL timing_start('iom_rstput')  
     392         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
    378393         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      ) 
    379          CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm  
     394         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
    380395         IF(nn_timing == 2)  CALL timing_stop('iom_rstput') 
    381396         ! 
    … …  
    385400         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    386401         CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 
    387          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )  
     402         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 
    388403      ENDIF 
    389404      !                       ! print mean trends (used for debugging) 
    390405      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
    391406      ! 
    392       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    393       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )  
     407      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
     408      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 
    394409      ! 
    395410      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr') 
    … …  
    407422      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio 
    408423      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The 
    409       !!      default values correspond to clear water (type I in Jerlov'  
     424      !!      default values correspond to clear water (type I in Jerlov' 
    410425      !!      (1968) classification. 
    411426      !!         called by tra_qsr at the first timestep (nit000) 
    … …  
    434449      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init') 
    435450      ! 
    436       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    437       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     451      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
     452      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
    438453      ! 
    439454 
    … …  
    465480 
    466481      IF( ln_traqsr ) THEN     ! control consistency 
    467          !                       
     482         ! 
    468483         IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio )   THEN 
    469484            CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 
    … …  
    480495            &              ' 2 bands, 3 RGB bands or bio-model light penetration' ) 
    481496         ! 
    482          IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1  
     497         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1 
    483498         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2 
    484499         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3 
    … …  
    498513      ENDIF 
    499514      !                          ! ===================================== ! 
    500       IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !   
     515      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  ! 
    501516         !                       ! ===================================== ! 
    502517         ! 
     518         ALLOCATE( xsi0r(jpi,jpj) ) 
    503519         xsi0r = 1.e0 / rn_si0 
    504520         xsi1r = 1.e0 / rn_si1 
    … …  
    539555                  zchl = 0.05                                 ! constant chlorophyll 
    540556                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    541                   zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll  
     557                  zekb(:,:) = rkrgb(1,irgb)                   ! Separation in R-G-B depending of the chlorophyll 
    542558                  zekg(:,:) = rkrgb(2,irgb) 
    543559                  zekr(:,:) = rkrgb(3,irgb) 
    … …  
    546562                  ze0(:,:,1) = rn_abs 
    547563                  ze1(:,:,1) = zcoef 
    548                   ze2(:,:,1) = zcoef  
     564                  ze2(:,:,1) = zcoef 
    549565                  ze3(:,:,1) = zcoef 
    550566                  zea(:,:,1) = tmask(:,:,1)                   ! = ( ze0+ze1+z2+ze3 ) * tmask 
    551                 
     567 
    552568                  DO jk = 2, nksr+1 
    553569!CDIR NOVERRCHK 
    554570                     DO jj = 1, jpj 
    555 !CDIR NOVERRCHK    
     571!CDIR NOVERRCHK 
    556572                        DO ji = 1, jpi 
    557                            zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r    ) 
     573                           zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) 
    558574                           zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 
    559575                           zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 
    … …  
    566582                        END DO 
    567583                     END DO 
    568                   END DO  
     584                  END DO 
    569585                  ! 
    570586                  DO jk = 1, nksr 
    … …  
    600616                  DO jj = 1, jpj                              ! top 400 meters 
    601617                     DO ji = 1, jpi 
    602                         zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    603                         zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    604                         etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)  
     618                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
     619                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
     620                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1) 
    605621                     END DO 
    606622                  END DO 
    … …  
    611627         ENDIF 
    612628         !                       ! ===================================== ! 
    613       ELSE                       !        No light penetration           !                    
     629      ELSE                       !        No light penetration           ! 
    614630         !                       ! ===================================== ! 
    615631         IF(lwp) THEN 
    … …  
    617633            WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 
    618634            WRITE(numout,*) '~~~~~~~~~~~~' 
    619             IF(lflush) CALL flush(numout) 
     635            IF(lflush) CALL flush(numout)             
    620636         ENDIF 
    621637      ENDIF 
    … …  
    630646      ENDIF 
    631647      ! 
    632       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    633       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     648      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
     649      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 
    634650      ! 
    635651      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_init') 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90

    r6486 r15603  
    2929   USE lib_mpp        ! MPP library 
    3030   USE wrk_nemo       ! Memory allocation 
     31   USE stopack      
    3132 
    3233   IMPLICIT NONE 
    … …  
    6869      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    6970      IF( ln_dyn_trd )   CALL trd_dyn_iom( putrd, pvtrd, ktrd, kt ) 
     71      IF( ln_dyn_trd .AND. ln_stopack .AND. ln_sppt_dyn ) & 
     72         & CALL dyn_sppt_collect( putrd, pvtrd, ktrd, kt ) 
    7073          
    7174      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r9163 r15603  
    3232   USE iom            ! I/O manager library 
    3333   USE lib_mpp        ! MPP library 
     34   USE stopack 
    3435   USE wrk_nemo       ! Memory allocation 
    3536 
    … …  
    271272      !                   ! 3D output of tracers trends using IOM interface 
    272273      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt ) 
    273  
    274       !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     274      IF( ln_tra_trd .AND. ln_stopack .AND. ln_sppt_tra )  & 
     275         & CALL tra_sppt_collect( ptrdx, ptrdy, ktrd, kt ) 
     276 
     277      !  Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    275278      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt ) 
    276279 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r12555 r15603  
    2626   USE wrk_nemo        ! Memory Allocation 
    2727   USE phycst, ONLY: vkarmn 
     28   USE stopack 
    2829 
    2930   IMPLICIT NONE 
    … …  
    5253   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
    5354   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction 
    54    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
     55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d,bfrcoef2d0   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
    5556 
    5657   !! * Substitutions 
    … …  
    6869      !!                ***  FUNCTION zdf_bfr_alloc  *** 
    6970      !!---------------------------------------------------------------------- 
    70       ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
     71      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), bfrcoef2d0(jpi,jpj),STAT=zdf_bfr_alloc ) 
    7172      ! 
    7273      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc ) 
    … …  
    107108         IF(lflush) CALL flush(numout) 
    108109      ENDIF 
     110      ! 
     111#if defined key_traldf_c2d || key_traldf_c3d 
     112      IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN 
     113         bfrcoef2d = bfrcoef2d0 
     114         CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr) 
     115      ENDIF 
     116#else 
     117      IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) & 
     118         & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// & 
     119                          'key_traldf_c2d or key_traldf_c3d') 
     120#endif 
    109121      ! 
    110122      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only 
    … …  
    135147               END DO 
    136148            END IF 
    137          !    
     149         ! 
    138150         ELSE 
    139151            zbfrt(:,:) = bfrcoef2d(:,:) 
    … …  
    178190               DO ji = 2, jpim1 
    179191                  ! (ISF) ======================================================================== 
    180                   ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     192                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
    181193                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points) 
    182194                  ! 
    … …  
    359371            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    360372         ENDIF 
    361           
     373 
    362374         IF ( ln_isfcav ) THEN 
    363375            IF(ln_tfr2d) THEN 
    … …  
    494506      ENDIF 
    495507      ! 
     508      bfrcoef2d0(:,:) = bfrcoef2d(:,:) 
    496509      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init') 
    497510      ! 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

    r12555 r15603  
    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 
    2728 
    2829   IMPLICIT NONE 
    … …  
    3031 
    3132   PUBLIC   zdf_evd    ! called by step.F90 
     33   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: rn_avevd0 
    3234 
    3335   !! * Substitutions 
    … …  
    4345      !!---------------------------------------------------------------------- 
    4446      !!                  ***  ROUTINE zdf_evd  *** 
    45       !!                    
     47      !! 
    4648      !! ** Purpose :   Local increased the vertical eddy viscosity and diffu- 
    4749      !!      sivity coefficients when a static instability is encountered. 
    4850      !! 
    4951      !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients 
    50       !!      are set to avevd (namelist parameter) if the water column is  
     52      !!      are set to avevd (namelist parameter) if the water column is 
    5153      !!      statically unstable (i.e. if rn2 < -1.e-12 ) 
    5254      !! 
    … …  
    7072         IF(lwp) WRITE(numout,*) 
    7173         IF(lwp .AND. lflush) CALL flush(numout) 
     74         ALLOCATE ( rn_avevd0(jpi,jpj) ) 
     75         rn_avevd0(:,:) = rn_avevd 
    7276      ENDIF 
    7377 
    7478      zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application 
     79 
     80#if defined key_traldf_c2d || key_traldf_c3d 
     81      IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN 
     82         rn_avevd0(:,:) = rn_avevd 
     83         CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd) 
     84      ENDIF 
     85#else 
     86      IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) & 
     87         & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// & 
     88                          'key_traldf_c2d or key_traldf_c3d') 
     89#endif 
    7590 
    7691      SELECT CASE ( nn_evdm ) 
    … …  
    8095         zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application 
    8196         ! 
    82          DO jk = 1, jpkm1  
     97         DO jk = 1, jpkm1 
    8398            DO jj = 2, jpj             ! no vector opt. 
    8499               DO ji = 2, jpi 
    … …  
    89104                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 
    90105#endif 
    91                      avt (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk) 
    92                      avm (ji  ,jj  ,jk) = rn_avevd * tmask(ji  ,jj  ,jk) 
    93                      avmu(ji  ,jj  ,jk) = rn_avevd * umask(ji  ,jj  ,jk) 
    94                      avmu(ji-1,jj  ,jk) = rn_avevd * umask(ji-1,jj  ,jk) 
    95                      avmv(ji  ,jj  ,jk) = rn_avevd * vmask(ji  ,jj  ,jk) 
    96                      avmv(ji  ,jj-1,jk) = rn_avevd * vmask(ji  ,jj-1,jk) 
     106                     avt (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk) 
     107                     avm (ji  ,jj  ,jk) = rn_avevd0(ji,jj) * tmask(ji  ,jj  ,jk) 
     108                     avmu(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji  ,jj  ,jk) 
     109                     avmu(ji-1,jj  ,jk) = rn_avevd0(ji,jj) * umask(ji-1,jj  ,jk) 
     110                     avmv(ji  ,jj  ,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj  ,jk) 
     111                     avmv(ji  ,jj-1,jk) = rn_avevd0(ji,jj) * vmask(ji  ,jj-1,jk) 
    97112                  ENDIF 
    98113               END DO 
    99114            END DO 
    100          END DO  
     115         END DO 
    101116         CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions 
    102117         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. ) 
    … …  
    105120         CALL iom_put( "avm_evd", zavm_evd )              ! output this change 
    106121         ! 
    107       CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12)  
     122      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
    108123         DO jk = 1, jpkm1 
    109 !!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL!  
     124!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
    110125            DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call) 
    111126               DO ji = 1, jpi 
    112127#if defined key_zdfkpp 
    113128                  ! no evd mixing in the boundary layer with KPP 
    114                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   &           
     129                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12  .AND.  fsdepw(ji,jj,jk) > hkpp(ji,jj)  )   & 
    115130#else 
    116131                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   & 
    117132#endif 
    118                      avt(ji,jj,jk) = rn_avevd * tmask(ji,jj,jk) 
     133                     avt(ji,jj,jk) = rn_avevd0(ji,jj) * tmask(ji,jj,jk) 
    119134               END DO 
    120135            END DO 
    121136         END DO 
    122137         ! 
    123       END SELECT  
     138      END SELECT 
    124139 
    125140      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r12555 r15603  
    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) 
     34   USE stopack 
    3435 
    3536   IMPLICIT NONE 
    … …  
    6162   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
    6263   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    63    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) 
    6465 
    6566   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
    66    REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1  
    67    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 
    6869   REAL(wp) ::   rghmin        = -0.28_wp 
    6970   REAL(wp) ::   rgh0          =  0.0329_wp 
    … …  
    7273   REAL(wp) ::   ra2           =  0.74_wp 
    7374   REAL(wp) ::   rb1           = 16.60_wp 
    74    REAL(wp) ::   rb2           = 10.10_wp          
    75    REAL(wp) ::   re2           =  1.33_wp          
     75   REAL(wp) ::   rb2           = 10.10_wp 
     76   REAL(wp) ::   re2           =  1.33_wp 
    7677   REAL(wp) ::   rl1           =  0.107_wp 
    7778   REAL(wp) ::   rl2           =  0.0032_wp 
    … …  
    133134      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments 
    134135      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars 
    135       REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      -  
     136      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
    136137      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      - 
    137138      REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      - 
    … …  
    139140      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep 
    140141      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar 
    141       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves  
     142      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
    142143      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves) 
    143144      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before 
    … …  
    156157      CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 
    157158      CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  ) 
    158        
     159 
    159160      ! Preliminary computing 
    160161 
    … …  
    165166         avm (:,:,:) = avm_k (:,:,:) 
    166167         avmu(:,:,:) = avmu_k(:,:,:) 
    167          avmv(:,:,:) = avmv_k(:,:,:)  
     168         avmv(:,:,:) = avmv_k(:,:,:) 
    168169      ENDIF 
    169170 
    170171      ! Compute surface and bottom friction at T-points 
    171 !CDIR NOVERRCHK           
    172       DO jj = 2, jpjm1           
    173 !CDIR NOVERRCHK          
    174          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. 
    175176            ! 
    176177            ! surface friction 
    177178            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
    178             !    
    179             ! bottom friction (explicit before friction)         
    180             ! Note that we chose here not to bound the friction as in dynbfr)    
    181             ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &          
    182                & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )       
    183             zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &          
    184                & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )       
    185             ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)          
    186          END DO          
    187       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 
    188189 
    189190      ! Set surface roughness length 
    190191      SELECT CASE ( nn_z0_met ) 
    191192      ! 
    192       CASE ( 0 )             ! Constant roughness           
     193      CASE ( 0 )             ! Constant roughness 
    193194         zhsro(:,:) = rn_hsro 
    194195      CASE ( 1 )             ! Standard Charnock formula 
    … …  
    226227      IF( nn_clos == 0 ) THEN    ! Mellor-Yamada 
    227228         DO jk = 2, jpkm1 
    228             DO jj = 2, jpjm1  
     229            DO jj = 2, jpjm1 
    229230               DO ji = fs_2, fs_jpim1   ! vector opt. 
    230231                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) 
    … …  
    275276               IF( ln_sigpsi ) THEN 
    276277                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1. 
    277                   zwall_psi(ji,jj,jk) = rsc_psi /   &  
     278                  zwall_psi(ji,jj,jk) = rsc_psi /   & 
    278279                     &     (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  ) 
    279280               ELSE 
    … …  
    294295               ! diagonal 
    295296               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  & 
    296                   &                       + rdt * zdiss * tmask(ji,jj,jk)  
     297                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
    297298               ! 
    298299               ! right hand side in en 
    … …  
    316317      ! First level 
    317318      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    318       en(:,:,1) = MAX(en(:,:,1), rn_emin)  
     319      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
    319320      z_elem_a(:,:,1) = en(:,:,1) 
    320321      z_elem_c(:,:,1) = 0._wp 
    321322      z_elem_b(:,:,1) = 1._wp 
    322       !  
     323      ! 
    323324      ! One level below 
    324325      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    325326          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    326327      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    327       z_elem_a(:,:,2) = 0._wp  
     328      z_elem_a(:,:,2) = 0._wp 
    328329      z_elem_c(:,:,2) = 0._wp 
    329330      z_elem_b(:,:,2) = 1._wp 
    … …  
    334335      ! Dirichlet conditions at k=1 
    335336      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    336       en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
     337      en(:,:,1)       = MAX(en(:,:,1), rn_emin) 
    337338      z_elem_a(:,:,1) = en(:,:,1) 
    338339      z_elem_c(:,:,1) = 0._wp 
    … …  
    357358      SELECT CASE ( nn_bc_bot ) 
    358359      ! 
    359       CASE ( 0 )             ! Dirichlet  
     360      CASE ( 0 )             ! Dirichlet 
    360361         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 
    361362         !                      ! Balance between the production and the dissipation terms 
    … …  
    377378               z_elem_c(ji,jj,ibotm1) = 0._wp 
    378379               z_elem_b(ji,jj,ibotm1) = 1._wp 
    379                en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin )  
     380               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
    380381            END DO 
    381382         END DO 
    382383         ! 
    383384      CASE ( 1 )             ! Neumman boundary condition 
    384          !                       
     385         ! 
    385386!CDIR NOVERRCHK 
    386387         DO jj = 2, jpjm1 
    … …  
    428429         END DO 
    429430      END DO 
    430       !                                            ! set the minimum value of tke  
     431      !                                            ! set the minimum value of tke 
    431432      en(:,:,:) = MAX( en(:,:,:), rn_emin ) 
    432433 
    … …  
    470471            DO jj = 2, jpjm1 
    471472               DO ji = fs_2, fs_jpim1   ! vector opt. 
    472                   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 
    473474               END DO 
    474475            END DO 
    … …  
    489490               ! 
    490491               ! psi / k 
    491                zratio = psi(ji,jj,jk) / eb(ji,jj,jk)  
     492               zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
    492493               ! 
    493494               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) 
    … …  
    509510               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term 
    510511               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 
    511                !                                                         
     512               ! 
    512513               ! building the matrix 
    513514               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 
    … …  
    551552      z_elem_c(:,:,2) = 0._wp 
    552553      z_elem_b(:,:,2) = 1._wp 
    553       !  
     554      ! 
    554555      ! 
    555556      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    … …  
    575576      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    576577 
    577       !    
     578      ! 
    578579      ! 
    579580      END SELECT 
    … …  
    585586      ! 
    586587      ! 
    587       CASE ( 0 )             ! Dirichlet  
     588      CASE ( 0 )             ! Dirichlet 
    588589         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 
    589590         !                      ! Balance between the production and the dissipation terms 
    … …  
    610611         ! 
    611612      CASE ( 1 )             ! Neumman boundary condition 
    612          !                       
     613         ! 
    613614!CDIR NOVERRCHK 
    614615         DO jj = 2, jpjm1 
    … …  
    692693            DO jj = 2, jpjm1 
    693694               DO ji = fs_2, fs_jpim1   ! vector opt. 
    694                   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) 
    695696               END DO 
    696697            END DO 
    … …  
    719720               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin ) 
    720721               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 
    721                ! 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) 
    722723               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    723724               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 
    724725            END DO 
    725726         END DO 
    726       END DO  
     727      END DO 
    727728 
    728729      ! 
    … …  
    806807            END DO 
    807808         END DO 
     809#if defined key_traldf_c2d || key_traldf_c3d 
     810         IF( ln_stopack) THEN 
     811            IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk) 
     812            IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk) 
     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 
    808819      END DO 
    809820      ! 
    … …  
    846857      !!---------------------------------------------------------------------- 
    847858      !!                  ***  ROUTINE zdf_gls_init  *** 
    848       !!                      
    849       !! ** Purpose :   Initialization of the vertical eddy diffivity and  
     859      !! 
     860      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
    850861      !!      viscosity when using a gls turbulent closure scheme 
    851862      !! 
    … …  
    881892      READ  ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 ) 
    882893902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp ) 
    883       IF(lwm .AND. nprint > 2) WRITE ( numond, namzdf_gls ) 
     894      IF(lwm) WRITE ( numond, namzdf_gls ) 
    884895 
    885896      IF(lwp) THEN                     !* Control print 
    … …  
    10691080      END SELECT 
    10701081      ! 
    1071       IF(lwp .AND. lflush) CALL flush(numout)  
     1082      IF(lwp .AND. lflush) CALL flush(numout) 
    10721083      !                                !* Set Schmidt number for psi diffusion in the wave breaking case 
    10731084      !                                     ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 
    10741085      !                                     !  or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 
    10751086      IF( ln_sigpsi ) THEN 
    1076          ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf  
     1087         ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 
    10771088         ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 
    10781089         ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 
    … …  
    10821093         rsc_psi0 = rsc_psi 
    10831094      ENDIF 
    1084   
     1095 
    10851096      !                                !* Shear free turbulence parameters 
    10861097      ! 
    … …  
    11291140      rc04  = rc03 * rc0 
    11301141      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1131       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
     1142      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
    11321143      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1133       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1144      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
    11341145      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    1135       rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
     1146      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 
    11361147      rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
    1137       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
     1148      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
    11381149 
    11391150      rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
    … …  
    11501161         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
    11511162      END DO 
    1152       !                               
     1163      ! 
    11531164      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files 
    11541165      ! 
    … …  
    11611172      !!--------------------------------------------------------------------- 
    11621173      !!                   ***  ROUTINE ts_rst  *** 
    1163       !!                      
     1174      !! 
    11641175      !! ** Purpose :   Read or write TKE file (en) in restart file 
    11651176      !! 
    11661177      !! ** Method  :   use of IOM library 
    1167       !!                if the restart does not contain TKE, en is either  
     1178      !!                if the restart does not contain TKE, en is either 
    11681179      !!                set to rn_emin or recomputed (nn_igls/=0) 
    11691180      !!---------------------------------------------------------------------- 
    … …  
    11771188      !!---------------------------------------------------------------------- 
    11781189      ! 
    1179       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     1190      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
    11801191         !                                   ! --------------- 
    11811192         IF( ln_rstart ) THEN                   !* Read the restart file 
    … …  
    11961207               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   ) 
    11971208               IF(nn_timing == 2)  CALL timing_stop('iom_rstget') 
    1198             ELSE                         
     1209            ELSE 
    11991210               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 
    12001211               IF(lwp .AND. lflush) CALL flush(numout) 
    12011212               en  (:,:,:) = rn_emin 
    1202                mxln(:,:,:) = 0.05         
     1213               mxln(:,:,:) = 0.05 
    12031214               avt_k (:,:,:) = avt (:,:,:) 
    12041215               avm_k (:,:,:) = avm (:,:,:) 
    … …  
    12111222            IF(lwp .AND. lflush) CALL flush(numout) 
    12121223            en  (:,:,:) = rn_emin 
    1213             mxln(:,:,:) = 0.05        
     1224            mxln(:,:,:) = 0.05 
    12141225         ENDIF 
    12151226         ! 
    … …  
    12191230         IF(lwp .AND. lflush) CALL flush(numout) 
    12201231         IF(nn_timing == 2)  CALL timing_start('iom_rstput') 
    1221          CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )  
     1232         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
    12221233         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  ) 
    12231234         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  ) 
    1224          CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )  
     1235         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
    12251236         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 
    12261237         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   ) 
    … …  
    12411252      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?' 
    12421253   END SUBROUTINE zdf_gls_init 
    1243     
     1254 
    12441255   SUBROUTINE zdf_gls( kt )          ! Empty routine 
    12451256   IMPLICIT NONE 
    1246       INTEGER, INTENT(in) ::   kt ! ocean time step    
     1257      INTEGER, INTENT(in) ::   kt ! ocean time step 
    12471258      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 
    12481259   END SUBROUTINE zdf_gls 
    1249     
     1260 
    12501261   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine 
    12511262   IMPLICIT NONE 
    … …  
    12581269   !!====================================================================== 
    12591270END MODULE zdfgls 
    1260  

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r12555 r15603  
    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 
    5757   USE agrif_opa_update 
    5858#endif 
    59  
    60  
     59   USE stopack 
    6160 
    6261   IMPLICIT NONE 
    … …  
    7574   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1) 
    7675   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
    77    REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation  
     76   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
    7877   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke 
    7978   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2] 
    … …  
    9493 
    9594   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    96    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
    97    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    9895   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9996#if defined key_c1d 
    … …  
    10299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers 
    103100#endif 
     101   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
     103   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    104104 
    105105   !! * Substitutions 
    … …  
    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  
     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 
     189      IF( ln_stopack) THEN 
     190         IF( nn_spp_tkelc > 0 ) THEN 
     191             rn_lc0 = rn_lc 
     192             CALL spp_gen( kt, rn_lc0, nn_spp_tkelc, rn_tkelc_sd, jk_spp_tkelc ) 
     193         ENDIF 
     194         IF( nn_spp_tkedf > 0 ) THEN 
     195             rn_ediff0 = rn_ediff 
     196             CALL spp_gen( kt, rn_ediff0, nn_spp_tkedf, rn_tkedf_sd, jk_spp_tkedf ) 
     197         ENDIF 
     198         IF( nn_spp_tkeds > 0 ) THEN 
     199             rn_ediss0 = rn_ediss 
     200             CALL spp_gen( kt, rn_ediss0, nn_spp_tkeds, rn_tkeds_sd, jk_spp_tkeds ) 
     201         ENDIF 
     202         IF( nn_spp_tkebb > 0 ) THEN 
     203             rn_ebb0 = rn_ebb 
     204             CALL spp_gen( kt, rn_ebb0, nn_spp_tkebb, rn_tkebb_sd, jk_spp_tkebb ) 
     205        ENDIF 
     206         IF( nn_spp_tkefr > 0 ) THEN 
     207             rn_efr0 = rn_efr 
     208             CALL spp_gen( kt, rn_efr0, nn_spp_tkefr, rn_tkefr_sd, jk_spp_tkefr ) 
     209         ENDIF 
     210      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 
    187216      ! 
    188217      CALL tke_tke      ! now tke (en) 
    … …  
    190219      CALL tke_avn      ! now avt, avm, avmu, avmv 
    191220      ! 
    192       avt_k (:,:,:) = avt (:,:,:)  
    193       avm_k (:,:,:) = avm (:,:,:)  
    194       avmu_k(:,:,:) = avmu(:,:,:)  
    195       avmv_k(:,:,:) = avmv(:,:,:)  
     221      avt_k (:,:,:) = avt (:,:,:) 
     222      avm_k (:,:,:) = avm (:,:,:) 
     223      avmu_k(:,:,:) = avmu(:,:,:) 
     224      avmv_k(:,:,:) = avmv(:,:,:) 
    196225      ! 
    197226#if defined key_agrif 
    198       ! Update child grid f => parent grid  
     227      ! Update child grid f => parent grid 
    199228      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    200 #endif       
    201      !  
     229#endif 
     230      IF (  kt == nitend ) THEN 
     231        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 
     232      ENDIF 
     233      ! 
     234 
    202235   END SUBROUTINE zdf_tke 
    203236 
    … …  
    225258      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
    226259      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
    227       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
    228       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
     260      REAL(wp) ::   zesh2                           ! temporary scalars 
     261      REAL(wp) ::   zfact1                          !    -         - 
    229262      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
    230263      REAL(wp) ::   ztau  , zdif                    !    -         - 
    … …  
    233266!!bfr      REAL(wp) ::   zebot                           !    -         - 
    234267      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
    235       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc 
     268      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3 
    236269      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 
    237270      !!-------------------------------------------------------------------- 
    … …  
    240273      ! 
    241274      CALL wrk_alloc( jpi,jpj, imlc )    ! integer 
    242       CALL wrk_alloc( jpi,jpj, zhlc )  
    243       CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    244       ! 
    245       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    246       zfact1 = -.5_wp * rdt  
    247       zfact2 = 1.5_wp * rdt * rn_ediss 
    248       zfact3 = 0.5_wp       * rn_ediss 
     275      CALL wrk_alloc( jpi,jpj, zhlc ) 
     276      CALL wrk_alloc( jpi,jpj, zbbrau ) 
     277      CALL wrk_alloc( jpi,jpj, zfact2 ) 
     278      CALL wrk_alloc( jpi,jpj, zfact3 ) 
     279      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
     280      ! 
     281      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation 
     282      zfact1 = -.5_wp * rdt 
     283      zfact2 = 1.5_wp * rdt * rn_ediss0 
     284      zfact3 = 0.5_wp       * rn_ediss0 
    249285      ! 
    250286      ! 
    … …  
    261297      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    262298         DO ji = fs_2, fs_jpim1   ! vector opt. 
    263             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    264          END DO 
    265       END DO 
    266        
     299            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 
     300         END DO 
     301      END DO 
     302 
    267303!!bfr   - start commented area 
    268304      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    303339         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    304340         DO jk = jpkm1, 2, -1 
    305             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 
    306342               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1) 
    307343                  zus  = zcof * taum(ji,jj) 
    … …  
    311347         END DO 
    312348         !                               ! finite LC depth 
    313          DO jj = 1, jpj  
     349         DO jj = 1, jpj 
    314350            DO ji = 1, jpi 
    315351               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 
    … …  
    326362                  !                                           ! vertical velocity due to LC 
    327363                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 
    328                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
     364                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    329365                  !                                           ! TKE Langmuir circulation source term 
    330                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   & 
     366                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) )* ( zwlc * zwlc * zwlc ) /   & 
    331367                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     368 
    332369               END DO 
    333370            END DO 
    … …  
    347384            DO ji = 1, jpi 
    348385               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    349                   &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
     386                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
    350387                  &                            / (  fse3uw_n(ji,jj,jk)               & 
    351388                  &                              *  fse3uw_b(ji,jj,jk)  ) 
    … …  
    376413                  !                                                           ! shear prod. at w-point weightened by mask 
    377414               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    378                   &    + ( 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) ) 
    379416                  ! 
    380417               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    381418               zd_lw(ji,jj,jk) = zzd_lw 
    382                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
     419               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
    383420               ! 
    384421               !                                   ! right hand side in en 
    385422               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    386                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     423                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    387424                  &                                 * wmask(ji,jj,jk) 
    388425            END DO 
    … …  
    433470      END DO 
    434471 
    435       !                                 ! Save TKE prior to nn_etau addition   
    436       e_niw(:,:,:) = en(:,:,:)   
    437       !   
     472      !                                 ! Save TKE prior to nn_etau addition 
     473      e_niw(:,:,:) = en(:,:,:) 
     474      ! 
    438475      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    439476      !                            !  TKE due to surface and internal wave breaking 
    … …  
    455492            DO jj = 2, jpjm1 
    456493               DO ji = fs_2, fs_jpim1   ! vector opt. 
    457                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     494                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    458495                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    459496               END DO 
    … …  
    464501            DO ji = fs_2, fs_jpim1   ! vector opt. 
    465502               jk = nmln(ji,jj) 
    466                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     503               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    467504                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    468505            END DO 
    … …  
    477514                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    478515                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    479                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    480                   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 
    481518                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    482                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     519                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    483520                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    484521               END DO 
    … …  
    486523         END DO 
    487524      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up) 
    488          IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau  
     525         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau 
    489526            DO jj = 2, jpjm1 
    490527               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    504541      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    505542      ! 
    506       DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term   
    507          DO jj = 2, jpjm1   
    508             DO ji = fs_2, fs_jpim1   ! vector opt.   
    509                e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)   
    510             END DO   
    511          END DO   
    512       END DO   
    513       !   
    514       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. ) 
    515552      ! 
    516553      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    517       CALL wrk_dealloc( jpi,jpj, zhlc )  
    518       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 ) 
    519559      ! 
    520560      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
    … …  
    529569      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity 
    530570      !! 
    531       !! ** Method  :   At this stage, en, the now TKE, is known (computed in  
    532       !!              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 
    533573      !!      computed from en and the strafification (N^2), then the mixings 
    534574      !!      coefficients are computed. 
    535575      !!              - Mixing length : a first evaluation of the mixing lengh 
    536576      !!      scales is: 
    537       !!                      mxl = sqrt(2*en) / N   
     577      !!                      mxl = sqrt(2*en) / N 
    538578      !!      where N is the brunt-vaisala frequency, with a minimum value set 
    539579      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean. 
    540       !!        The mixing and dissipative length scale are bound as follow :  
     580      !!        The mixing and dissipative length scale are bound as follow : 
    541581      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom. 
    542582      !!                        zmxld = zmxlm = mxl 
    543583      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 
    544       !!         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 
    545585      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 
    546586      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings 
    547       !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to  
    548       !!                    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 
    549589      !!                    scales are: 
    550       !!                        zmxld = sqrt( lup * ldown )  
     590      !!                        zmxld = sqrt( lup * ldown ) 
    551591      !!                        zmxlm = min ( lup , ldown ) 
    552592      !!              - Vertical eddy viscosity and diffusivity: 
    553593      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 
    554       !!                      avt = max( avmb, pdlr * avm )   
     594      !!                      avt = max( avmb, pdlr * avm ) 
    555595      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 
    556596      !! 
    … …  
    567607      IF( nn_timing == 1 )  CALL timing_start('tke_avn') 
    568608 
    569       CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
     609      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
    570610 
    571611      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    … …  
    576616      ! 
    577617      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
    578       zmxlm(:,:,:)  = rmxl_min     
     618      zmxlm(:,:,:)  = rmxl_min 
    579619      zmxld(:,:,:)  = rmxl_min 
    580620      ! 
    … …  
    586626            END DO 
    587627         END DO 
    588       ELSE  
     628      ELSE 
    589629         zmxlm(:,:,1) = rn_mxl0 
    590630      ENDIF 
    … …  
    604644      !                     !* Physical limits for the mixing length 
    605645      ! 
    606       zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value  
     646      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
    607647      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    608648      ! 
    … …  
    698738            DO ji = fs_2, fs_jpim1   ! vector opt. 
    699739               zsqen = SQRT( en(ji,jj,jk) ) 
    700                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     740               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 
    701741               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    702742               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    … …  
    704744            END DO 
    705745         END DO 
     746#if defined key_traldf_c2d || key_traldf_c3d 
     747         IF( ln_stopack) THEN 
     748            IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 
     749            IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 
     750         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 
    706756      END DO 
    707757      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    … …  
    749799      ENDIF 
    750800      ! 
    751       CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )  
     801      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
    752802      ! 
    753803      IF( nn_timing == 1 )  CALL timing_stop('tke_avn') 
    … …  
    759809      !!---------------------------------------------------------------------- 
    760810      !!                  ***  ROUTINE zdf_tke_init  *** 
    761       !!                      
    762       !! ** Purpose :   Initialization of the vertical eddy diffivity and  
     811      !! 
     812      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
    763813      !!              viscosity when using a tke turbulent closure scheme 
    764814      !! 
    … …  
    776826         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   & 
    777827         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   & 
    778          &                 nn_etau , nn_htau  , rn_efr , rn_c    
     828         &                 nn_etau , nn_htau  , rn_efr , rn_c 
    779829      !!---------------------------------------------------------------------- 
    780830 
    … …  
    843893         rn_mxl0 = rmxl_min 
    844894      ENDIF 
    845        
     895 
     896      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc 
     897      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 
     898      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 
     899      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb 
     900      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr 
     901 
    846902      IF( nn_etau == 2  ) THEN 
    847903          ierr = zdf_mxl_alloc() 
    … …  
    855911 
    856912      !                               !* depth of penetration of surface tke 
    857       IF( nn_etau /= 0 ) THEN       
     913      IF( nn_etau /= 0 ) THEN 
    858914         htau(:,:) = 0._wp 
    859915         SELECT CASE( nn_htau )             ! Choice of the depth of penetration 
    … …  
    861917            htau(:,:) = 10._wp 
    862918         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees 
    863             htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )             
     919            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   ) 
    864920         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer 
    865921            rhtau = -1._wp / LOG( 1._wp - rn_c ) 
    … …  
    904960      END DO 
    905961      dissl(:,:,:) = 1.e-12_wp 
    906       !                               
     962      ! 
    907963      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files 
    908964      ! 
    … …  
    913969     !!--------------------------------------------------------------------- 
    914970     !!                   ***  ROUTINE tke_rst  *** 
    915      !!                      
     971     !! 
    916972     !! ** Purpose :   Read or write TKE file (en) in restart file 
    917973     !! 
    918974     !! ** Method  :   use of IOM library 
    919      !!                if the restart does not contain TKE, en is either  
    920      !!                set to rn_emin or recomputed  
     975     !!                if the restart does not contain TKE, en is either 
     976     !!                set to rn_emin or recomputed 
    921977     !!---------------------------------------------------------------------- 
    922978     INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    … …  
    927983     !!---------------------------------------------------------------------- 
    928984     ! 
    929      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     985     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
    930986        !                                   ! --------------- 
    931987        IF( ln_rstart ) THEN                   !* Read the restart file 
    … …  
    9561012              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
    9571013              ! 
    958               avt_k (:,:,:) = avt (:,:,:) 
    959               avm_k (:,:,:) = avm (:,:,:) 
    960               avmu_k(:,:,:) = avmu(:,:,:) 
    961               avmv_k(:,:,:) = avmv(:,:,:) 
    962               ! 
    9631014              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    9641015           ENDIF 
    9651016        ELSE                                   !* Start from rest 
    9661017           en(:,:,:) = rn_emin * tmask(:,:,:) 
    967            DO jk = 1, jpk                           ! set the Kz to the background value 
    968               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    969               avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    970               avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
    971               avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    972            END DO 
    9731018        ENDIF 
    9741019        ! 
     1020        avt_k (:,:,:) = avt (:,:,:) 
     1021        avm_k (:,:,:) = avm (:,:,:) 
     1022        avmu_k(:,:,:) = avmu(:,:,:) 
     1023        avmv_k(:,:,:) = avmv(:,:,:) 
     1024 
    9751025     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    9761026        !                                   ! ------------------- 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r13070 r15603  
    8585   USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges  
    8686   USE sbc_oce, ONLY: lk_oasis 
     87   USE stopack 
    8788   USE stopar 
    8889   USE stopts 
    … …  
    498499                            CALL dia_hsb_init   ! heat content, salt content and volume budgets 
    499500                            CALL     trd_init   ! Mixed-layer/Vorticity/Integral constraints trends 
    500                             CALL     bias_init  ! Pressure correction bias 
     501                             CALL stopack_init   ! STOPACK scheme                             
     502                            CALL    bias_init   ! Pressure correction bias 
    501503      IF( lk_diaobs     ) THEN                  ! Observation & model comparison 
    502504                            CALL dia_obs_init            ! Initialize observational data 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/step.F90

    r14987 r15603  
    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 
    … …  
    119119      END DO 
    120120 
     121      IF( ln_stopack )   CALL stopack_pert( kstp ) 
    121122                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    122123                                                      ! clem: moved here for bdy ice purpose 
     124 
    123125      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    124126      ! Update stochastic parameters and random T/S fluctuations 
    125127      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    126        IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters 
    127        IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations 
     128 
     129      IF( ln_sto_eos )                 CALL sto_par( kstp )   ! Stochastic parameters 
     130      IF( ln_sto_eos )                 CALL sto_pts( tsn  )   ! Random T/S fluctuations 
     131      IF( ln_stopack .AND. ln_skeb  )  CALL skeb_comp( kstp ) ! Stochastic Energy Backscatter 
    128132 
    129133      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    149153      ENDIF 
    150154      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
    151          DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
     155#if defined key_traldf_c2d || key_traldf_c3d 
     156         IF ( ln_stopack .AND. ( nn_spp_arnf > 0 ) ) THEN 
     157              rn_avt_rnf0 = rn_avt_rnf 
     158              CALL spp_gen( kstp, rn_avt_rnf0,nn_spp_arnf,rn_arnf_sd,jk_spp_arnf ) 
     159         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 
     165         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf0(:,:) * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO 
    152166      ENDIF 
    153167      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity 
    … …  
    163177      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' ) 
    164178      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' ) 
     179      IF( lrst_oce .AND. ln_stopack)   CALL stopack_rst( kstp, 'WRITE' ) 
    165180      ! 
    166181      !  LATERAL  PHYSICS 
    … …  
    195210      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    196211                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur) 
    197       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    198                          CALL wzv           ( kstp )  ! now cross-level velocity  
    199  
    200       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 
    201216          ! In case the time splitting case, update almost all momentum trends here: 
    202217          ! Note that the computation of vertical velocity above, hence "after" sea level 
    203218          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
     219            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations 
    204220                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    205221            IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
    … …  
    214230                                  va(:,:,:) = 0.e0 
    215231          IF(  lk_asminc .AND. ln_asmiau .AND. & 
    216              & ln_dyninc       )  CALL dyn_asm_inc  ( kstp )   ! apply dynamics assimilation increment 
    217           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! subtract Neptune velocities (simplified) 
    218           IF( lk_bdy           )  CALL bdy_dyn3d_dmp( kstp )   ! bdy damping trends 
    219                                   CALL dyn_adv      ( kstp )   ! advection (vector or flux form) 
    220                                   CALL dyn_vor      ( kstp )   ! vorticity term including Coriolis 
    221                                   CALL dyn_ldf      ( kstp )   ! lateral mixing 
    222           IF( ln_neptsimp )       CALL dyn_nept_cor ( kstp )   ! add Neptune velocities (simplified) 
     232             & ln_dyninc       )  CALL dyn_asm_inc   ( kstp )   ! apply dynamics assimilation increment 
     233          IF( ln_stopack .AND. ln_sppt_dyn ) & 
     234             &                    CALL dyn_sppt_apply( kstp, 0 ) 
     235          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! subtract Neptune velocities (simplified) 
     236          IF( lk_bdy           )  CALL bdy_dyn3d_dmp ( kstp )   ! bdy damping trends 
     237                                  CALL dyn_adv       ( kstp )   ! advection (vector or flux form) 
     238                                  CALL dyn_vor       ( kstp )   ! vorticity term including Coriolis 
     239                                  CALL dyn_ldf       ( kstp )   ! lateral mixing 
     240          IF( ln_neptsimp )       CALL dyn_nept_cor  ( kstp )   ! add Neptune velocities (simplified) 
    223241#if defined key_agrif 
    224242          IF(.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn        ! momentum sponge 
    … …  
    232250                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 
    233251          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    234                                   CALL wzv           ( kstp )  ! now cross-level velocity  
     252                                  CALL wzv           ( kstp )  ! now cross-level velocity 
    235253      ENDIF 
    236254 
    … …  
    262280                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero 
    263281 
    264       IF(  lk_asminc .AND. ln_asmiau .AND. & 
    265          & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment 
     282      IF( ln_stopack .AND. ln_sppt_tra ) & 
     283         &                   CALL tra_sppt_apply ( kstp, 0 ) 
     284      IF( lk_asminc .AND. ln_asmiau .AND. & 
     285        & ln_trainc      )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment 
    266286                             CALL tra_sbc    ( kstp )       ! surface boundary condition 
    267287      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
    … …  
    280300      IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge 
    281301#endif 
     302      IF( ln_stopack .AND. ln_sppt_tra ) & 
     303         &                   CALL tra_sppt_apply ( kstp, 1 ) 
    282304                             CALL tra_zdf    ( kstp )       ! vertical mixing and after tracer fields 
    283305 
    … …  
    285307         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
    286308                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
     309         IF( ln_stopack .AND. ln_sppt_tra ) & 
     310            &                CALL tra_sppt_apply ( kstp, 2 ) 
    287311                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    288312            IF( ln_zps .AND. .NOT. ln_isfcav)                                & 
    … …  
    300324               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    301325               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    302          IF( ln_zps .AND.       ln_isfcav)                                   &  
     326         IF( ln_zps .AND.       ln_isfcav)                                   & 
    303327               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    304328               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    … …  
    308332                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    309333         IF( ln_bias )       CALL dyn_bias( kstp ) 
     334         IF( ln_stopack .AND. ln_sppt_tra )   THEN 
     335                                CALL tra_sppt_apply ( kstp, 2 ) 
     336                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
     337                                IF( ln_zps .AND. .NOT. ln_isfcav)                & 
     338                                & CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     339                                &                   rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     340                                IF( ln_zps .AND.       ln_isfcav)                                   & 
     341                                & CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     342                                &                   rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     343                                &                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
     344         ENDIF 
    310345      ENDIF 
    311346 
    … …  
    318353                               ua(:,:,:) = ua_sv(:,:,:) 
    319354                               va(:,:,:) = va_sv(:,:,:) 
    320                                                              ! Revert now divergence and rotational to previously computed ones  
     355                                                             ! Revert now divergence and rotational to previously computed ones 
    321356                                                             !(needed because of the time swap in div_cur, at the beginning of each time step) 
    322357                               hdivn(:,:,:) = hdivb(:,:,:) 
    323                                rotn(:,:,:)  = rotb(:,:,:)  
     358                               rotn(:,:,:)  = rotb(:,:,:) 
    324359 
    325360                               CALL dyn_bfr( kstp )         ! bottom friction 
     361            IF( ln_stopack .AND. ln_sppt_dyn ) & 
     362               &               CALL dyn_sppt_apply ( kstp, 1 ) 
    326363                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    327364      ELSE 
    … …  
    331368        IF(  lk_asminc .AND. ln_asmiau .AND. & 
    332369           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
     370        IF( ln_stopack .AND. ln_sppt_dyn ) & 
     371           &                   CALL dyn_sppt_apply ( kstp, 0 ) 
    333372        IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
    334373        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
    … …  
    343382                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    344383                               CALL dyn_bfr( kstp )         ! bottom friction 
     384            IF( ln_stopack .AND. ln_sppt_dyn ) & 
     385               &               CALL dyn_sppt_apply ( kstp, 1 ) 
    345386                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    346387                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    347388      ENDIF 
    348389                               CALL dyn_nxt( kstp )         ! lateral velocity at next time step 
     390            IF( ln_stopack .AND. ln_sppt_dyn ) & 
     391               &               CALL dyn_sppt_apply ( kstp, 2 ) 
     392            IF( ln_stopack .AND. ln_skeb ) & 
     393               &               CALL skeb_apply ( kstp ) 
    349394 
    350395                               CALL ssh_swp( kstp )         ! swap of sea surface height 
    … …  
    352397      ! 
    353398 
    354       ! Moved bias _wrt to before rst_write as used restart parameters for nn_stocklist option that are changed by rst_wrt    
     399      ! Moved bias write to before rst_write 
    355400      IF( lrst_bias )          CALL bias_wrt     ( kstp ) 
    356  
     401       
    357402      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
    358403      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
    … …  
    361406      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    362407      ! AGRIF 
    363       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    364                                CALL Agrif_Integrate_ChildGrids( stp )   
     408      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     409                               CALL Agrif_Integrate_ChildGrids( stp ) 
    365410 
    366411      IF ( Agrif_NbStepint().EQ.0 ) THEN 
    … …  
    393438      ! 
    394439#if defined key_iomput 
    395       IF( kstp == nitend .OR. indic < 0 ) THEN  
     440      IF( kstp == nitend .OR. indic < 0 ) THEN 
    396441                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    397          IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !  
     442         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
    398443      ENDIF 
    399444#endif 
    400445      ! 
    401446      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset 
    402       !      
     447      ! 
    403448      ! 
    404449   END SUBROUTINE stp 

  • branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r8400 r15603  
    9696   USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine) 
    9797   USE diaharm 
     98   USE diaprod          ! ocean model: product diagnostics 
    9899   USE flo_oce          ! floats variables 
    99100   USE floats           ! floats computation               (flo_stp routine) 
    … …  
    113114   USE timing           ! Timing 
    114115 
     116   USE stopack          ! Stochastic physics 
     117 
    115118#if defined key_agrif 
    116119   USE agrif_opa_sponge ! Momemtum and tracers sponges 
Note: See TracChangeset for help on using the changeset viewer.