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 11334 for NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90 – NEMO

Ignore:
Timestamp:
2019-07-24T11:20:42+02:00 (5 years 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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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 
Note: See TracChangeset for help on using the changeset viewer.