Changeset 11334


Ignore:
Timestamp:
2019-07-24T11:20:42+02:00 (13 months ago)
Author:
flemarie
Message:

Follow-on to the ABL implementation (see ticket #2131)

  • Modification of diawri.F90 to output ABL variables with IOIPSL (i.e. without key_iomput)
  • Modification of ice_var_agg in icevar.F90 to allow the use of tm_su in ABL and bulk
  • Error handling in case ln_abl = .true. and ABL sources were not compiled
  • ABL now works with SI3 considering an average over ice categories
  • Update reference namelist to avoid runtime warnings due to nam_tide
  • Sanity checks with ORCA2 for the following configurations :

1 - ABL src + ln_blk = .true.
2 - ABL src + ln_abl = .true.
3 - no ABL src + ln_blk = .true.

All configurations are MPP reproducible and configurations 1 and 3 results are identical

Location:
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref

    r11322 r11334  
    301301   rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented) 
    302302   rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities) 
    303    ln_smth_pblh  = .false.    !  Smoothing of PBL height with a 2d Hanning filter 
     303   ln_smth_pblh  = .true.    !  Smoothing of PBL height with a 2d Hanning filter 
    304304/ 
    305305!----------------------------------------------------------------------- 
     
    589589!----------------------------------------------------------------------- 
    590590   ln_tide     = .false.      ! Activate tides 
    591       ln_tide_pot   = .true.                !  use tidal potential forcing 
    592          ln_scal_load  = .false.               ! Use scalar approximation for 
    593             rn_scal_load = 0.094               !     load potential 
    594          ln_read_load  = .false.               ! Or read load potential from file 
    595             cn_tide_load = 'tide_LOAD_grid_T.nc'  ! filename for load potential 
    596             !       
    597       ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup 
    598          rdttideramp   =    0.                 !  ramp duration in days 
    599       clname(1)     = 'DUMMY'               !  name of constituent - all tidal components must be set in namelist_cfg 
     591   ln_tide_pot   = .true.                !  use tidal potential forcing 
     592   ln_scal_load  = .false.               ! Use scalar approximation for  
     593   ln_read_load  = .false.               ! Or read load potential from file 
     594   cn_tide_load  = 'tide_LOAD_grid_T.nc'  ! filename for load potential      
     595   ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup 
     596   rn_scal_load  = 0.094               !     load potential  
     597   rdttideramp   =    0.                 !  ramp duration in days 
     598   clname        = 'FOUR'               !  name of constituent - all tidal components must be set in namelist_cfg 
    600599/ 
    601600!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90

    r11322 r11334  
    3838!=================================================================================================== 
    3939   SUBROUTINE abl_stp( kt, psst, pssu, pssv, pssq, &                            ! in 
    40               &            pu_dta, pv_dta, pt_dta, pq_dta, & 
    41               &            pslp_dta, pgu_dta, pgv_dta, & 
    42               &            pcd_du, psen, pevp,     &                            ! in/out 
    43               &            pwndm, ptaui, ptauj, ptaum ) 
     40              &            pu_dta, pv_dta, pt_dta, pq_dta,    &  
     41              &            pslp_dta, pgu_dta, pgv_dta,        & 
     42              &            pcd_du, psen, pevp,                &                            ! in/out 
     43              &            pwndm, ptaui, ptauj, ptaum         &  
     44#if defined key_si3           
     45              &          , ptm_su,pssu_ice,pssv_ice,pssq_ice,pcd_du_ice  &  
     46              &          , psen_ice, pevp_ice, pwndm_ice, pfrac_oce      &  
     47           &          , ptaui_ice,     ptauj_ice                      & 
     48#endif            
     49           &   ) 
    4450!--------------------------------------------------------------------------------------------------- 
    4551 
     
    8086      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaum      ! ||tau||  
    8187      ! 
     88#if defined key_si3     
     89      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptm_su       ! ice-surface temperature [K] 
     90      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssu_ice     ! ice-surface u (U-point) 
     91      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssv_ice     ! ice-surface v (V-point)       
     92      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pssq_ice     ! ice-surface humidity  
     93      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pcd_du_ice   ! Cd x Du over ice (T-point) 
     94     REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   psen_ice     ! Ch x Du over ice (T-point) 
     95     REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point) 
     96     REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||   
     97     REAL(wp) , INTENT(inout), DIMENSION(:,:  ) ::   pfrac_oce    ! 
     98      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point) 
     99      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)      
     100#endif      
     101     ! 
    82102      REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zrhoa, zwnd_i, zwnd_j 
    83 !      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   zFC 
    84103      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF     
    85 !      REAL(wp), DIMENSION(1:jpi,  jptq  )        ::   zBC 
    86       REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)     ::   z_cft      !--FL--to be removed after the test phase    
     104      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase    
    87105      ! 
    88106      REAL(wp), DIMENSION(1:jpi,1:jpka  )        ::   z_elem_a 
     
    91109      ! 
    92110      INTEGER             ::   ji, jj, jk, jtra, jbak               ! dummy loop indices 
    93       REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1 
    94       REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv 
     111      REAL(wp)            ::   zztmp, zcff, ztemp, zhumi, zcff1, zztmp1, zztmp2 
     112      REAL(wp)            ::   zcff2, zfcor, zmsk, zsig, zcffu, zcffv, zzice,zzoce 
    95113      ! 
    96114      !!---------------------------------------------------------------------       
     
    108126      DO jj = 1,jpj 
    109127         DO ji = 1,jpi 
    110             ustar2(ji,jj) = pCd_du(ji,jj) * pwndm(ji,jj) 
     128          zzoce         = pCd_du    (ji,jj) * pwndm    (ji,jj) 
     129#if defined key_si3            
     130         zzice         = pCd_du_ice(ji,jj) * pwndm_ice(ji,jj)  
     131          ustar2(ji,jj) = zzoce * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * zzice  
     132#else 
     133           ustar2(ji,jj) = zzoce  
     134#endif                   
    111135         END DO 
    112136      END DO   
     
    154178            IF(jtra == jp_ta) THEN 
    155179               DO ji = 1,jpi  ! surface boundary condition for temperature               
    156                   z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * psen(ji, jj)  
    157                   tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   & 
    158                     &                               + rdt * psen(ji, jj) * ( psst(ji, jj) + rt0 )                 
    159                   tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
     180                  zztmp1 = psen(ji, jj)  
     181              zztmp2 = psen(ji, jj) * ( psst(ji, jj) + rt0 )    
     182#if defined key_si3              
     183                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj)  
     184                  zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * psen_ice(ji,jj) * ptm_su(ji,jj) 
     185#endif               
     186              z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1             
     187                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2                
     188              tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    160189               END DO  
    161190            ELSE    
    162191               DO ji = 1,jpi  ! surface boundary condition for humidity               
    163                   z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * pevp(ji, jj)  
    164                   tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra )   & 
    165                     &                               + rdt * pevp(ji, jj) * pssq(ji, jj)                
     192                  zztmp1 = pevp(ji, jj)  
     193              zztmp2 = pevp(ji, jj) * pssq(ji, jj)    
     194#if defined key_si3              
     195                  zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji,jj)  
     196              zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pevp_ice(ji, jj) * pssq_ice(ji, jj)     
     197#endif    
     198                  z_elem_b( ji,     2              ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1 
     199                  tq_abl  ( ji, jj, 2, nt_a, jtra  ) = e3t_abl( 2 ) * tq_abl  ( ji, jj, 2, nt_n, jtra ) + rdt * zztmp2                
    166200                  tq_abl  ( ji, jj, jpka, nt_a, jtra ) = e3t_abl( jpka ) * tq_abl  ( ji, jj, jpka, nt_n, jtra ) 
    167201               END DO                                      
     
    251285         ENDDO  
    252286      END IF 
    253  
     287      
    254288      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    255289      !                            !  4 *** Advance u,v to time n+1  
     
    273307            z_elem_a( ji,     2    ) = 0._wp 
    274308            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                        
    275             zcff  = rdt * pcd_du(ji, jj) 
    276             z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff  
    277             u_abl( ji, jj,    2, nt_a ) = u_abl( ji, jj,    2, nt_a ) + 0.5_wp * zcff * ( pssu(ji-1, jj) + pssu(ji,jj) )  
     309            ! 
     310         zztmp1  = pcd_du(ji, jj) 
     311         zztmp2  = 0.5_wp * pcd_du(ji, jj) * ( pssu(ji-1, jj) + pssu(ji,jj) )        
     312#if defined key_si3              
     313            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 
     314         zzice  = 0.5_wp * ( pssu_ice(ji-1, jj) + pssu_ice(ji,jj) ) 
     315         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
     316#endif            
     317         z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1           
     318            u_abl( ji, jj,    2, nt_a ) =      u_abl( ji, jj,    2, nt_a ) + rdt * zztmp2 
     319          
    278320            !++ Top Neumann B.C. 
    279321            !z_elem_a( ji,     jpka ) = - 0.5_wp * rdt * ( Avm_abl( ji, jj, jpka )+ Avm_abl( ji+1, jj, jpka ) ) / e3w_abl( jpka )  
     
    331373            !++ Surface boundary condition 
    332374            z_elem_a( ji,     2    ) = 0._wp 
    333             z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )                                        
    334             zcff  = rdt * pcd_du(ji, jj) 
    335             z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + zcff  
    336             v_abl( ji, jj,    2, nt_a ) =  v_abl( ji, jj, 2, nt_a ) + 0.5_wp * zcff * ( pssv(ji, jj) + pssv(ji, jj-1) )   
     375            z_elem_c( ji,     2    ) = - rdt * Avm_abl( ji, jj, 2   ) / e3w_abl( 2   )         
     376            ! 
     377         zztmp1 = pcd_du(ji, jj) 
     378         zztmp2 = 0.5_wp * pcd_du(ji, jj) * ( pssv(ji, jj) + pssv(ji, jj-1) )   
     379#if defined key_si3              
     380            zztmp1 = zztmp1 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) 
     381         zzice  = 0.5_wp * ( pssv_ice(ji, jj) + pssv_ice(ji,jj-1) ) 
     382         zztmp2 = zztmp2 * pfrac_oce(ji,jj) + (1._wp - pfrac_oce(ji,jj)) * pcd_du_ice(ji, jj) * zzice 
     383#endif          
     384            z_elem_b( ji,     2       ) = e3t_abl( 2 ) - z_elem_c( ji, 2 ) + rdt * zztmp1  
     385            v_abl( ji, jj,    2, nt_a ) =         v_abl( ji, jj, 2, nt_a ) + rdt * zztmp2 
    337386            !++ Top Neumann B.C.             
    338387            !z_elem_a( ji,     jpka ) = -rdt * Avm_abl( ji, jj, jpka ) / e3w_abl( jpka )  
     
    372421      END DO             ! end outer loop  
    373422      !-------------    
    374    
     423     
    375424      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    376425      !                            !  5 *** Apply nudging on the dynamics and the tracers  
     
    428477      !------------- 
    429478      END DO             ! end outer loop 
    430       !-------------  
     479      !-------------       
    431480      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    432481      !                            !  6 *** MPI exchanges   
    433482      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    434       ! 
    435       CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), kfillmode = jpfillnothing  ) 
     483     ! 
     484      CALL lbc_lnk_multi( 'ablmod', u_abl(:,:,:,nt_a), 'T', -1., v_abl(:,:,:,nt_a), 'T', -1., kfillmode = jpfillnothing  ) 
    436485      CALL lbc_lnk_multi( 'ablmod', tq_abl(:,:,:,nt_a,jp_ta), 'T', 1., tq_abl(:,:,:,nt_a,jp_qa), 'T', 1., kfillmode = jpfillnothing  )   ! ++++ this should not be needed... 
    437       ! 
     486     ! 
    438487      CALL iom_put ( 'u_abl',  u_abl(:,:,2:jpka,nt_a      )) 
    439488      CALL iom_put ( 'v_abl',  v_abl(:,:,2:jpka,nt_a      )) 
     
    504553         CALL prt_ctl( tab2d_2=ptauj  , clinfo2=          'vtau   : ' ) 
    505554      ENDIF 
    506       
     555 
     556#if defined key_si3 
     557         ! ------------------------------------------------------------ ! 
     558         !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     559         ! ------------------------------------------------------------ ! 
     560         DO jj = 2, jpjm1 
     561            DO ji = 2, jpim1   
     562                
     563            zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     564            zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     565             
     566            ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     567                  &                      +    zrhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     568                  &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
     569               ptauj_ice(ji,jj) = 0.5_wp * (  zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
     570                  &                      +    zrhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
     571                  &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
     572            END DO 
     573         END DO 
     574         CALL lbc_lnk_multi( 'ablmod', ptaui_ice, 'U', -1., ptauj_ice, 'V', -1. ) 
     575         ! 
     576         IF(ln_ctl)   CALL prt_ctl( tab2d_1=ptaui_ice  , clinfo1=' abl_stp: putaui : '   & 
     577            &                     , tab2d_2=ptauj_ice  , clinfo2='          pvtaui : ' ) 
     578#endif 
    507579      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    508580      !                            !  8 *** Swap time indices for the next timestep 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90

    r11322 r11334  
    3030   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3131   USE prtctl         ! Print control 
    32  
     32#if defined key_si3 
     33   USE ice    , ONLY : u_ice, v_ice, tm_su, ato_i      ! ato_i = total open water fractional area 
     34   USE sbc_ice, ONLY : wndm_ice, utau_ice, vtau_ice 
     35#endif    
     36#if ! defined key_iomput 
     37   USE diawri    , ONLY : dia_wri_alloc_abl 
     38#endif 
    3339   IMPLICIT NONE 
    3440   PRIVATE 
     
    141147      CALL iom_get( inum, jpdom_unknown, 'ghw_abl', ghw_abl(:) ) 
    142148      CALL iom_close( inum ) 
     149      
     150#if ! defined key_iomput 
     151     IF( dia_wri_alloc_abl()  /= 0 ) CALL ctl_stop( 'STOP', 'abl_init : unable to allocate arrays' ) 
     152#endif 
    143153 
    144154      IF(lwp) THEN 
     
    303313      !! 
    304314      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp 
     315#if defined key_si3 
     316      REAL(wp), DIMENSION(jpi,jpj) ::   zssqi, zcd_dui, zseni, zevpi 
     317#endif      
    305318      INTEGER                      ::   jbak, jbak_dta, ji, jj 
    306319      !!--------------------------------------------------------------------- 
     
    319332              &           sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & 
    320333              &              zssq, zcd_du, zsen, zevp ) 
    321    
     334 
     335#if defined key_si3 
     336     CALL blk_ice_1( u_abl(:,:,2,nt_n), v_abl(:,:,2,nt_n), tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   & 
     337        &            sf(jp_slp)%fnow(:,:,1), u_ice, v_ice,                                          & 
     338         &            pseni=zseni, pevpi=zevpi, ptsui=tm_su, pssqi=zssqi, pcd_dui=zcd_dui )   ! outputs   
     339#endif   
     340 
    322341      !!------------------------------------------------------------------------------------------- 
    323342      !! 3 - Advance ABL variables from now (n) to after (n+1) 
     
    329348              &         sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),                          & 
    330349              &         zcd_du, zsen, zevp,    &                                ! in/out 
    331               &         wndm, utau, vtau, taum )                                ! out 
    332  
     350              &         wndm, utau, vtau, taum &                                 ! out 
     351#if defined key_si3           
     352              &          , tm_su, u_ice, v_ice, zssqi, zcd_dui   &  
     353              &          , zseni, zevpi, wndm_ice, ato_i,        &  
     354           &            utau_ice,     vtau_ice                & 
     355#endif            
     356              &                                                                ) 
    333357      !!------------------------------------------------------------------------------------------- 
    334358      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since  
     
    340364              &           sst_m , zsen, zevp ) 
    341365 
     366#if defined key_si3 
     367! Avoid a USE abl in icesbc module 
     368      sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa) 
     369#endif  
     370 
    342371   END SUBROUTINE sbc_abl 
    343372 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icesbc.F90

    r11322 r11334  
    7575            &                                      u_ice, v_ice,   &    ! inputs 
    7676            &                                      putaui = utau_ice, pvtaui = vtau_ice )    ! outputs                              
    77          CASE( jp_abl     )   ;    CALL blk_ice_1( sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   & 
    78             &                                      sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1),   & 
    79             &                                      u_ice, v_ice,   &    ! inputs 
    80             &                                      putaui = utau_ice, pvtaui = vtau_ice )    ! outputs   
     77 !        CASE( jp_abl     )    utau_ice & vtau_ice are computed in ablmod 
    8178         CASE( jp_purecpl )   ;    CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )   ! Coupled      formulation 
    8279      END SELECT 
     
    146143      CASE( jp_usr )              !--- user defined formulation 
    147144                                  CALL usrdef_sbc_ice_flx( kt, h_s, h_i ) 
    148       CASE( jp_blk, jp_abl )  !--- bulk formulation (temporary version for ABL option) 
     145      CASE( jp_blk, jp_abl )  !--- bulk formulation & ABL formulation 
    149146                                  CALL blk_ice_2    ( t_su, h_s, h_i, alb_ice, sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),    & 
    150147            &                                           sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) )    !  
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90

    r11223 r11334  
    114114      ! 
    115115      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction   
    116  
     116      ! 
     117      ALLOCATE( z1_at_i(jpi,jpj) ) 
     118      WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
     119      ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
     120      END WHERE 
     121      ! 
     122      tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
     123      WHERE( at_i(:,:)<=epsi20 ); tm_su(:,:) = rt0; END WHERE       
     124      ! 
    117125      ! The following fields are calculated for diagnostics and outputs only 
    118126      ! ==> Do not use them for other purposes 
    119127      IF( kn > 1 ) THEN 
    120128         ! 
    121          ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
    122          WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
    123          ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
    124          END WHERE 
     129         ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
    125130         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:) 
    126131         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp 
     
    135140         !          
    136141         !                          ! mean temperature (K), salinity and age 
    137          tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    138142         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    139143         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
     
    153157         !                           ! put rt0 where there is no ice 
    154158         WHERE( at_i(:,:)<=epsi20 ) 
    155             tm_su(:,:) = rt0 
    156159            tm_si(:,:) = rt0 
    157160            tm_i (:,:) = rt0 
     
    159162         END WHERE 
    160163 
    161          DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 
     164         DEALLOCATE( z1_vt_i , z1_vt_s ) 
    162165      ENDIF 
     166      ! 
     167      DEALLOCATE( z1_at_i ) 
    163168      ! 
    164169   END SUBROUTINE ice_var_agg 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90

    r10425 r11334  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
     28   USE abl            ! abl variables in case ln_abl = .true. 
    2829   USE dom_oce        ! ocean space and time domain 
    2930   USE phycst         ! physical constants 
     
    6667   PUBLIC   dia_wri_state 
    6768   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
    68  
     69#if ! defined key_iomput    
     70   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.) 
     71#endif 
    6972   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    70    INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
     73   INTEGER ::          nb_T, ndim_bT                 ! grid_T file 
    7174   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file 
    7275   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7376   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
     77   INTEGER ::   ndim_A, ndim_hA   ! grid_T file    
     78   INTEGER ::   nid_A, nz_A, nh_A                    ! grid_ABL file    
    7479   INTEGER ::   ndex(1)                              ! ??? 
    75    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
    76    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
     80   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV, ndex_hA 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V, ndex_A 
    7782   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    7883 
     
    409414         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    410415         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
    411          ! 
    412       dia_wri_alloc = MAXVAL(ierr) 
     416     dia_wri_alloc = MAXVAL(ierr) 
    413417      CALL mpp_sum( 'diawri', dia_wri_alloc ) 
    414418      ! 
    415419   END FUNCTION dia_wri_alloc 
    416  
    417     
     420  
     421   INTEGER FUNCTION dia_wri_alloc_abl() 
     422      !!---------------------------------------------------------------------- 
     423     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 
     424      CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 
     425      ! 
     426   END FUNCTION dia_wri_alloc_abl 
     427  
    418428   SUBROUTINE dia_wri( kt ) 
    419429      !!--------------------------------------------------------------------- 
     
    435445      INTEGER  ::   ji, jj, jk                               ! dummy loop indices 
    436446      INTEGER  ::   ierr                                     ! error code return from allocation 
    437       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers 
     447      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma, ipka   ! local integers 
    438448      INTEGER  ::   jn, ierror                               ! local integers 
    439449      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars 
     
    441451      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    442452      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     453      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! 3D workspace 
    443454      !!---------------------------------------------------------------------- 
    444455      !  
     
    449460         ninist = 0 
    450461      ENDIF 
     462      
    451463      ! 
    452464      ! 0. Initialisation 
     
    472484      ijmi = 1      ;      ijma = jpj 
    473485      ipk = jpk 
     486     IF(ln_abl) ipka = jpkam1 
    474487 
    475488      ! define time axis 
     
    574587            &          "m", ipk, gdepw_1d, nz_W, "down" ) 
    575588 
     589         IF( ln_abl ) THEN  
     590         ! Define the ABL grid FILE ( nid_A ) 
     591            CALL dia_nam( clhstnam, nwrite, 'grid_ABL' ) 
     592            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     593            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     594               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     595               &          nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 
     596            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept 
     597               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 
     598            !                                                            ! Index of ocean points 
     599         ALLOCATE( zw3d_abl(jpi,jpj,ipka) )  
     600         zw3d_abl(:,:,:) = 1._wp  
     601         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume 
     602            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface 
     603         DEALLOCATE(zw3d_abl) 
     604         ENDIF 
    576605 
    577606         ! Declare all the output fields as NETCDF variables 
     
    623652         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm 
    624653            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     654! 
     655         IF( ln_abl ) THEN 
     656         !                                                                                      !!! nid_A : 3D 
     657         CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
     658               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     659            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl 
     660               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     661            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl 
     662               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     663            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl 
     664               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     665            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl 
     666               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     667            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl 
     668               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     669            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl 
     670               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     671            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh 
     672               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                  
     673#if defined key_si3 
     674            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i 
     675               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout ) 
     676#endif 
     677          CALL histend( nid_A, snc4chunks=snc4set ) 
     678       ! 
     679       ENDIF 
    625680! 
    626681         IF( ln_icebergs ) THEN 
     
    676731          
    677732         clmx ="l_max(only(x))"    ! max index on a period 
    678 !         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
    679 !            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout ) 
     733         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
     734            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout ) 
    680735#if defined key_diahth 
    681736         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth 
     
    790845      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction    
    791846      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed    
     847! 
     848      IF( ln_abl ) THEN  
     849        ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
     850        IF( ln_mskland )   THEN  
     851          DO jk=1,jpka 
     852             zw3d_abl(:,:,jk) = tmask(:,:,1) 
     853            END DO 
     854       ELSE 
     855            zw3d_abl(:,:,:) = 1._wp      
     856         ENDIF        
     857       CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
     858        CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
     859        CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
     860        CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
     861        CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl      
     862        CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
     863        CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
     864        CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl   
     865#if defined key_si3 
     866         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i 
     867#endif 
     868       DEALLOCATE(zw3d_abl) 
     869     ENDIF 
    792870! 
    793871      IF( ln_icebergs ) THEN 
     
    826904         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    827905      ENDIF 
    828 !      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
    829 !      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
     906      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
     907      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
    830908 
    831909#if defined key_diahth 
     
    862940         CALL histclo( nid_V ) 
    863941         CALL histclo( nid_W ) 
     942         IF(ln_abl) CALL histclo( nid_A ) 
    864943      ENDIF 
    865944      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcabl.F90

    r11275 r11334  
    88   !!---------------------------------------------------------------------- 
    99   USE sbc_oce, ONLY :   ght_abl, ghw_abl, e3t_abl, e3w_abl 
     10   USE lib_mpp, ONLY :   ctl_stop 
    1011 
    1112   IMPLICIT NONE 
     
    3132      !! 
    3233      !!---------------------------------------------------------------------- 
    33        
     34      CALL ctl_stop( 'STOP', 'ln_abl = .true. but ABL source directory was not included',   & 
     35        & '(Either switch to ln_abl = .false. or modify your cfg.txt file and recompile)' ) 
     36     !! 
    3437   END SUBROUTINE sbc_abl_init 
    3538 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90

    r11322 r11334  
    816816      ! 
    817817      INTEGER  ::   ji, jj    ! dummy loop indices 
    818       REAL(wp) ::   zwndi_t , zwndj_t            ! relative wind components at T-point 
     818      REAL(wp) ::   zwndi_t , zwndj_t, zootm_su, zztmp1, zztmp2 ! relative wind components at T-point 
    819819      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    820820      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui     ! transfer coefficient for momentum      (tau) 
     
    855855      zcd_dui(:,:) = wndm_ice(:,:) * Cd_atm(:,:) 
    856856       
    857  !     IF( ln_blk ) THEN  
     857      IF( ln_blk ) THEN  
    858858         ! ------------------------------------------------------------ ! 
    859859         !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     
    874874         IF(ln_ctl)   CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
    875875            &                     , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
    876  !     ELSE 
    877  !        DO jj = 1, jpj 
    878  !           DO ji = 1, jpi 
    879  !              pcd_dui(ji,jj) = zcd_dui(ji,jj) 
    880  !              pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_atm(ji,jj) 
    881  !              pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_atm(ji,jj) 
    882  !              pssqi  (ji,jj) = 11637800.0_wp * EXP( -5897.8_wp / (ptsui(ji,jj)+rt0) ) / zrhoa(ji,jj) 
    883  !           END DO 
    884  !        END DO 
    885  !     ENDIF 
     876      ELSE 
     877         zztmp1 = 11637800.0_wp 
     878       zztmp2 =    -5897.8_wp 
     879       DO jj = 1, jpj 
     880            DO ji = 1, jpi 
     881               pcd_dui(ji,jj) = zcd_dui (ji,jj) 
     882               pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_atm(ji,jj) 
     883               pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_atm(ji,jj) 
     884            zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??) 
     885            pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj) 
     886            END DO 
     887         END DO 
     888      ENDIF 
    886889      ! 
    887890      IF(ln_ctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
     
    12481251      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
    12491252      REAL(wp) ::   zCdn_form_tmp 
    1250       !!---------------------------------------------------------------------- 
    1251  
     1253      !!--------------------------------------------------------------------------- 
     1254      ! fl: global variable tm_su could be used directly instead of recomputing it  
     1255     ! 
    12521256      ! mean temperature 
    12531257      WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 
Note: See TracChangeset for help on using the changeset viewer.