












      subroutine bulk_flux (tile)      


      use Agrif_Util
      interface
        subroutine Sub_Loop_bulk_flux(tile,N2d,A2d,Mmmpi,Lmmpi,N3d)
          implicit none
      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      integer(4) :: N2d
      integer(4) :: Mmmpi
      integer(4) :: Lmmpi
      integer(4) :: N3d
      real, dimension(1:N2d,1:NSA,0:NPP-1) :: A2d
      integer(4) :: tile
        end subroutine Sub_Loop_bulk_flux

      end interface
      integer(4) :: tile
      

        call Sub_Loop_bulk_flux(tile, Agrif_tabvars_i(15)%iarray0, Agrif
     &_tabvars(117)%array3, Agrif_tabvars_i(8)%iarray0, Agrif_tabvars_i(
     &7)%iarray0, Agrif_tabvars_i(16)%iarray0)

      end


      subroutine Sub_Loop_bulk_flux(tile,N2d,A2d,Mmmpi,Lmmpi,N3d)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      integer(4) :: N2d
      integer(4) :: Mmmpi
      integer(4) :: Lmmpi
      integer(4) :: N3d
      real, dimension(1:N2d,1:NSA,0:NPP-1) :: A2d
      integer(4) :: tile

                   

                                              
      integer(4) :: trd
      integer(4) :: omp_get_thread_num
                                            
                             
                                            

                                                             
                                                            
                                           
                                                               
                                    
                                      
                                                          

                       

                                    

                        
                              

                                          
                                             
                          

                        

                                          
                        

                                                        
                      
                              

                                       
                                                                     
                                          
                         

                             

                             

                            

                            

                                             

                            

                                     

                                                               

                         

                                                      
                                                      
                                                                  
                               

                               

                              

                              

                               

                              

                              

                                

                             

                                                    
                                 
                                      
                                      
                                                            
      integer(4) :: chunk_size_X
      integer(4) :: margin_X
      integer(4) :: chunk_size_E
      integer(4) :: margin_E
                                             
      integer(4) :: Istr
      integer(4) :: Iend
      integer(4) :: Jstr
      integer(4) :: Jend
      integer(4) :: i_X
      integer(4) :: j_E
      chunk_size_X=(Lmmpi+NSUB_X-1)/NSUB_X
      margin_X=(NSUB_X*chunk_size_X-Lmmpi)/2
      chunk_size_E=(Mmmpi+NSUB_E-1)/NSUB_E
      margin_E=(NSUB_E*chunk_size_E-Mmmpi)/2
      j_E=tile/NSUB_X
      i_X=tile-j_E*NSUB_X
      Istr=1+i_X*chunk_size_X-margin_X
      Iend=Istr+chunk_size_X-1
      Istr=max(Istr,1)
      Iend=min(Iend,Lmmpi)
      Jstr=1+j_E*chunk_size_E-margin_E
      Jend=Jstr+chunk_size_E-1
      Jstr=max(Jstr,1)
      Jend=min(Jend,Mmmpi)
      trd=omp_get_thread_num()
      call bulk_flux_tile (Istr,Iend,Jstr,Jend,
     &                               A2d(1,1,trd),A2d(1,2,trd))
      return
       
      

      end subroutine Sub_Loop_bulk_flux

      subroutine bulk_flux_tile (Istr,Iend,Jstr,Jend, aer,cer)      


      use Agrif_Util
      interface
        subroutine Sub_Loop_bulk_flux_tile(Istr,Iend,Jstr,Jend,aer,cer,p
     &add_E,Mm,padd_X,Lm,vmask,vwnd,svstr,umask,uwnd,sustr,shflx_rlw,shf
     &lx_sen,shflx_lat,shflx_rsw,rmask,prate,srflx,stflx,radlw,rhum,nrhs
     &,t,tair,wspd,rho0,NORTH_INTER,SOUTH_INTER,EAST_INTER,WEST_INTER)
          implicit none
      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      real, parameter :: qp2 = 0.0000172D0
      real, parameter :: pi = 3.14159265358979323846D0
      real, parameter :: deg2rad = pi/180.D0
      real, parameter :: rad2deg = 180.D0/pi
      real, parameter :: Eradius = 6371315.0D0
      real, parameter :: Erotation = 7.292115090D-5
      real, parameter :: day2sec = 86400.D0
      real, parameter :: sec2day = 1.D0/86400.D0
      real, parameter :: year2day = 365.25D0
      real, parameter :: day2year = 1.D0/365.25D0
      real, parameter :: jul_off = 2440000.D0
      real, parameter :: g = 9.81D0
      real, parameter :: Cp = 3985.0D0
      real, parameter :: vonKar = 0.41D0
      real, parameter :: spval = -999.0D0
      logical, parameter :: mask_val = .true.
      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      INTEGER(4), parameter :: IterFl = 3
      INTEGER(4), parameter :: mb = 1
      REAL, parameter :: blk_ZW = 10.0D0
      REAL, parameter :: blk_ZT = 2.0D0
      REAL, parameter :: blk_ZToZW = blk_ZT/blk_ZW
      REAL, parameter :: Log10oLogZw = LOG(10.0D0*10000.D0)/LOG(blk_ZW*1
     &0000.D0)
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4) :: nrhs
      real :: rho0
      logical :: NORTH_INTER
      logical :: SOUTH_INTER
      logical :: EAST_INTER
      logical :: WEST_INTER
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: vmask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: vwnd
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: svstr
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: umask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: uwnd
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: sustr
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_rlw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_sen
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_lat
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_rsw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: rmask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: prate
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: srflx
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E,1:NT) :: stflx
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: radlw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: rhum
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E,1:N,1:3,1:NT) :: t
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: tair
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: wspd
      INTEGER(4) :: Istr
      INTEGER(4) :: Iend
      INTEGER(4) :: Jstr
      INTEGER(4) :: Jend
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: aer
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: cer
        end subroutine Sub_Loop_bulk_flux_tile

      end interface
      INTEGER(4) :: Istr
      INTEGER(4) :: Iend
      INTEGER(4) :: Jstr
      INTEGER(4) :: Jend
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: aer
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: cer
      

        call Sub_Loop_bulk_flux_tile(Istr,Iend,Jstr,Jend,aer,cer, Agrif_
     &tabvars_i(14)%iarray0, Agrif_tabvars_i(4)%iarray0, Agrif_tabvars_i
     &(13)%iarray0, Agrif_tabvars_i(2)%iarray0, Agrif_tabvars(78)%array2
     &, Agrif_tabvars(232)%array2, Agrif_tabvars(264)%array2, Agrif_tabv
     &ars(79)%array2, Agrif_tabvars(233)%array2, Agrif_tabvars(265)%arra
     &y2, Agrif_tabvars(248)%array2, Agrif_tabvars(246)%array2, Agrif_ta
     &bvars(247)%array2, Agrif_tabvars(249)%array2, Agrif_tabvars(81)%ar
     &ray2, Agrif_tabvars(237)%array2, Agrif_tabvars(215)%array2, Agrif_
     &tabvars(250)%array3, Agrif_tabvars(236)%array2, Agrif_tabvars(238)
     &%array2, Agrif_tabvars_i(154)%iarray0, Agrif_tabvars(148)%array5, 
     &Agrif_tabvars(239)%array2, Agrif_tabvars(234)%array2, Agrif_tabvar
     &s_r(12)%array0, Agrif_tabvars_l(13)%larray0, Agrif_tabvars_l(14)%l
     &array0, Agrif_tabvars_l(11)%larray0, Agrif_tabvars_l(12)%larray0)

      end


      subroutine Sub_Loop_bulk_flux_tile(Istr,Iend,Jstr,Jend,aer,cer,pad
     &d_E,Mm,padd_X,Lm,vmask,vwnd,svstr,umask,uwnd,sustr,shflx_rlw,shflx
     &_sen,shflx_lat,shflx_rsw,rmask,prate,srflx,stflx,radlw,rhum,nrhs,t
     &,tair,wspd,rho0,NORTH_INTER,SOUTH_INTER,EAST_INTER,WEST_INTER)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      integer(4), parameter :: LLm0 = 41
      integer(4), parameter :: MMm0 = 42
      integer(4), parameter :: N = 32
      integer(4), parameter :: NP_XI = 1
      integer(4), parameter :: NP_ETA = 4
      integer(4), parameter :: NNODES = NP_XI*NP_ETA
      integer(4), parameter :: NPP = 1
      integer(4), parameter :: NSUB_X = 1
      integer(4), parameter :: NSUB_E = 1
      integer(4), parameter :: NWEIGHT = 1000
      integer(4), parameter :: stdout = 6
      integer(4), parameter :: Np = N+1
      integer(4), parameter :: NSA = 28
      real, parameter :: Vtransform = 2
      integer(4), parameter :: itemp = 1
      integer(4), parameter :: ntrc_temp = 1
      integer(4), parameter :: ntrc_salt = 1
      integer(4), parameter :: ntrc_pas = 0
      integer(4), parameter :: ntrc_bio = 0
      integer(4), parameter :: ntrc_subs = 0
      integer(4), parameter :: ntrc_substot = 0
      integer(4), parameter :: ntrc_sed = 0
      integer(4), parameter :: NTA = itemp+ntrc_salt
      integer(4), parameter :: NT = itemp+ntrc_salt+ntrc_pas+ntrc_bio+nt
     &rc_sed
      integer(4), parameter :: NTot = NT
      integer(4), parameter :: isalt = itemp+1
      integer(4), parameter :: ntrc_diabio = 0
      integer(4), parameter :: ntrc_diats = 0
      integer(4), parameter :: ntrc_diauv = 0
      integer(4), parameter :: ntrc_diavrt = 0
      integer(4), parameter :: ntrc_diaek = 0
      integer(4), parameter :: ntrc_diapv = 0
      integer(4), parameter :: ntrc_diaeddy = 0
      integer(4), parameter :: ntrc_surf = 0
      real, parameter :: qp2 = 0.0000172D0
      real, parameter :: pi = 3.14159265358979323846D0
      real, parameter :: deg2rad = pi/180.D0
      real, parameter :: rad2deg = 180.D0/pi
      real, parameter :: Eradius = 6371315.0D0
      real, parameter :: Erotation = 7.292115090D-5
      real, parameter :: day2sec = 86400.D0
      real, parameter :: sec2day = 1.D0/86400.D0
      real, parameter :: year2day = 365.25D0
      real, parameter :: day2year = 1.D0/365.25D0
      real, parameter :: jul_off = 2440000.D0
      real, parameter :: g = 9.81D0
      real, parameter :: Cp = 3985.0D0
      real, parameter :: vonKar = 0.41D0
      real, parameter :: spval = -999.0D0
      logical, parameter :: mask_val = .true.
      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      INTEGER(4), parameter :: IterFl = 3
      INTEGER(4), parameter :: mb = 1
      REAL, parameter :: blk_ZW = 10.0D0
      REAL, parameter :: blk_ZT = 2.0D0
      REAL, parameter :: blk_ZToZW = blk_ZT/blk_ZW
      REAL, parameter :: Log10oLogZw = LOG(10.0D0*10000.D0)/LOG(blk_ZW*1
     &0000.D0)
      integer(4) :: padd_E
      integer(4) :: Mm
      integer(4) :: padd_X
      integer(4) :: Lm
      integer(4) :: nrhs
      real :: rho0
      logical :: NORTH_INTER
      logical :: SOUTH_INTER
      logical :: EAST_INTER
      logical :: WEST_INTER
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: vmask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: vwnd
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: svstr
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: umask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: uwnd
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: sustr
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_rlw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_sen
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_lat
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: shflx_rsw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: rmask
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: prate
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: srflx
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E,1:NT) :: stflx
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: radlw
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: rhum
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E,1:N,1:3,1:NT) :: t
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: tair
      real, dimension(-1:Lm+2+padd_X,-1:Mm+2+padd_E) :: wspd
      INTEGER(4) :: Istr
      INTEGER(4) :: Iend
      INTEGER(4) :: Jstr
      INTEGER(4) :: Jend
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: aer
      REAL, dimension(Istr-2:Iend+2,Jstr-2:Jend+2) :: cer

                   

                                            
                             
                                            

                                                             
                                                            
                                           
                                                               
                                    
                                      
                                                          

                       

                                    

                        
                              

                                          
                                             
                          

                        

                                          
                        

                                                        
                      
                              

                                       
                                                                     
                                          
                         

                             

                             

                            

                            

                                             

                            

                                     

                                                               

                         

                                                      
                                                      
                                                                  
                               

                               

                              

                              

                               

                              

                              

                                

                             

                                            
                                               
                                            
                                               
                                                                 
                                                 
                                 
                                               
                                               
                                               
                                               
                                               
                                               
                                             
                                             
                                             
                                             
                                             
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                               
                                              
                                                  
                                                  
                                                  
                                                  
                                                 
                                                 
                                               
                                               
                                                      
                                                 
                                                 
                                                 
                                                 
                                                 
                                                 
                                                 
                                                           
                                                           
                                                           
                                    
                                                
                                                
                                                
                                                
                                                 
                           
                           
                           
                           
                             
                                              
                          
                                                
                                                
                                                   
                                              
                                               
                                                   
                                                
                                                  
                                                 
                                                 
                                                              
                             
                                                 
                                                  
                                                 
                                                
                                             
                                                
                            
               
                                 

                                                
                                                
                                                     
                                                   
                                                   
                                                         
                                                
                                   
                                                    
                                              
                                     
                                                
                                              
                                            
                                                           
                                                     
                           
                                  
                                                
                                                
                                                     
                                                   
                                                   
                                                         
                                                           
                                                         
                                                              
                                                           
                                                                
                                                                   
                                                                    
                                                                   
                                                   
                                 
                                                    
                                      
                                                    
                                      
                                                    
                                      
                                                    
                                      
                                                      
                                   
                                        
                                        
                                                       
                                                     
                                             
                                                               
                                              
                                                           
                                        
                                                   
                                 
                                               
                                               
                                                
                                                
                                                
                                               
                                               
                                               
                              
                              
                                
                                
                                
                              
                              
                              
                                                  
                                                  
                                                   
                                                   
                                                   
                                                  
                                                  
                                                  
                                  
                                  
                                    
                                    
                                    
                               
                               
                                  
                                                              
                                
                                       
                                                           
                                                                
                                             
                                                  
                                                     
                           
                                                                       
                                                                       
                                                                 
                             
                                                                    
                                       
                                                              
                                                             
                                                                       
                                                                 
                                                      
                                        
                                                
                                 
                                                   
                                   
                                 
                                
                                           
                                         
                                             
                                                              
                                              
                                                                       
                                                                  
                                                                       
                                                                 
                                 
                                
                                                                       
                                                                 
                                                                 
                                                                 
                                                        
                                                                       
                                                                 
                                                                 
            
                                              
                                                   
                     
                              
                               
                                 
                                                  
                                                                       
                                                                 
                                                                 
                        
                      
                           
                                                                       
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                                 
                                                    
                           
                                      
                             
                                                 
                                                                       
                                                                 
                     
                                                      
                                                                       
                                                                 
              
                                          
                                                                       
                                                                 
       
                                                           
      real(8) :: Cu_Nbq_X
      real(8) :: Cu_Nbq_Y
      real(8) :: Cu_Nbq_Z
                                             
                                                                       
                                                                 
              
                                                  
                                                                       
                                                                 
              
                                  
                                           
                                                  
                                                                   
                                                               
                                                      
                                                                       
                                                                 
                      
                                                                       
                                                                 
                                                                 
                                                                 
                                                                 
                                                        
                                
                                                                       
                                                                 
                                                        

                                                                       
                                                                 
                  
                                                                       
                                                                 
                                                                 
                                                                 
                  

                          

              
                             

                  
                               

                 
                                

                       
                                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                              
      INTEGER(4) :: i
      INTEGER(4) :: j
      INTEGER(4) :: m
      INTEGER(4) :: iter
                                                       
                                                       
      INTEGER(4) :: imin
      INTEGER(4) :: imax
      INTEGER(4) :: jmin
      INTEGER(4) :: jmax
                                           
                                     
                                    
      REAL :: cff
      REAL :: cff2
                                                             
      REAL :: rho0i
      REAL :: cpi
      REAL :: RH
      REAL :: Hlv
      REAL :: wspd0
                                                              
      REAL :: TseaC
      REAL :: TseaK
      REAL :: Qsea
      REAL :: spec_hum
      REAL :: Tvstar
                                                          
      REAL :: TairC
      REAL :: TairK
      REAL :: Q
      REAL :: rhoAir
      REAL :: Bf
                                                                       
                                                                 
  
      REAL, dimension(1:2) :: delW
      REAL :: delT
      REAL :: delQ
      REAL, dimension(1:2) :: ddelW
      REAL :: ddelT
      REAL :: ddelQ
                                                                     
      REAL :: Cd
      REAL :: Ch
      REAL :: Ce
      REAL :: qsat
      REAL :: iexns
      REAL :: iexna
                                                                       
                                                                 
  
      REAL, dimension(1:2) :: Wstar
      REAL :: Tstar
      REAL :: Qstar
      REAL :: Wstar0
      REAL :: Tstar0
      REAL :: Qstar0
                                                                  
      REAL, dimension(1:2) :: psi_u
      REAL, dimension(1:2) :: ZoLu
      REAL :: psi_t
      REAL :: ZoLt
      REAL :: patm
                                                            
      REAL :: hfsen
      REAL :: hflat
      REAL :: hflw
      REAL :: upvel
      REAL :: evap
                                           
                                           
                                                    
                                          
      REAL :: logus10
      REAL :: logts10
                                                                       
                                                                 
                                                  
                                                                   
      REAL :: bulk_psiu_coare
      REAL :: bulk_psit_coare
      REAL :: air_visc
                                                               
      REAL :: iZoW
      REAL :: iZoT
      REAL :: CC
      REAL :: Ch10
      REAL :: Ribcu
      REAL :: charn
                                                        
      REAL :: iZo10
      REAL :: iZoT10
      REAL :: Ri
      REAL :: Rr
      REAL :: VisAir
                                                          
                                                          
      if (.not.WEST_INTER) then
        imin=Istr-1
      else
        imin=Istr-2
      endif
      if (.not.EAST_INTER) then
        imax=Iend+1
      else
        imax=Iend+2
      endif
      if (.not.SOUTH_INTER) then
        jmin=Jstr-1
      else
        jmin=Jstr-2
      endif
      if (.not.NORTH_INTER) then
        jmax=Jend+1
      else
        jmax=Jend+2
      endif
      rho0i=1.0D0/rho0
      cpi=1.0D0/cp
      iexns  = (psurf*ip00)**(-rdocpd)
      DO j=jmin,jmax
        DO i=imin,imax
          wspd0     = wspd(i,j)
          wspd0     = MAX ( wspd0 , 0.1D0 * MIN(10.D0, blk_ZW) )
          TairC = tair(i,j)
          TairK = TairC + CtoK
          TseaC = t(i,j,N,nrhs,itemp)
          TseaK = TseaC + CtoK
          RH     = rhum(i,j)
          Q      = spec_hum (RH,psurf,TairC)
          CALL exner_patm_from_tairabs(iexna,patm,Q,TairK,blk_ZT,psurf)
          rhoAir = patm*(1.D0+Q) / ( blk_Rgas*TairK*(1.D0+MvoMa*Q) )
          Qsea   = qsat(TseaK,psurf,0.98D0)
          delW(1) = SQRT(wspd0*wspd0+0.25D0)
          delQ = Q-Qsea
          cff  = CtoK*(iexna-iexns)
          delT = TairC*iexna - TseaC*iexns + cff
          Wstar(1:mb)  = 0.035D0*delW(1:mb)*Log10oLogZw
          VisAir = air_visc ( TairC )
          charn  = 0.011D0
          Ch10   = 0.00115D0
          Ribcu  = -    blk_ZW / ( blk_Zabl * 0.004D0 * blk_beta**3 )
          DO m=1,mb
            iZo10  = g*Wstar(m) /
     &               (charn*Wstar(m)*Wstar(m)*Wstar(m)+0.11D0*g*VisAir)
            iZoT10 = 0.1D0 * exp(vonKar*vonKar /
     &                     ( Ch10*LOG( 10.0D0*iZo10 ) ) )
            CC     = LOG( blk_ZW*iZo10 )*LOG( blk_ZW*iZo10 )/
     &               LOG( blk_ZT*iZoT10 )
            Ri     =  g * blk_ZW * ( delT+cpvir*TairK*delQ )/
     &                             ( TairK*delW(m)*delW(m) )
            IF ( Ri < 0.0D0 ) THEN
              ZoLu(m)=CC*Ri/(1.0D0+Ri/Ribcu)
            ELSE
              ZoLu(m)=CC*Ri/(1.0D0+3.0D0*Ri/CC)
            ENDIF
            psi_u(m)     = bulk_psiu_coare ( ZoLu(m) )
            logus10   = LOG(blk_ZW* iZo10)
            Wstar(m)  = delW(m)*vonKar/(logus10-psi_u(m))
          ENDDO
          ZoLt      = ZoLu(mb)*blk_ZToZW
          psi_t     = bulk_psit_coare ( ZoLt    )
          logts10   = LOG(blk_ZT*iZoT10)
          cff       =      vonKar/(logts10-psi_t)
          Tstar     = delT*cff
          Qstar     = delQ*cff
          charn = 0.011D0
          IF     ( delW(1) > 18.0D0 ) then
            charn = 0.018D0
          ELSEIF ( delW(1) > 10.0D0 ) then
            charn = 0.011D0+0.125D0*(0.018D0-0.011D0)*(delW(1)-10.D0)
          ENDIF
          DO Iter=1,IterFl
            DO m=1,mb
              iZoW    = g*Wstar(m) / ( charn*Wstar(m)*Wstar(m)*Wstar(m)
     &                                 +0.11D0*g*VisAir )
              Rr      = Wstar(m)/(iZow*VisAir)
              iZoT    = MAX(8695.65D0,18181.8D0*(Rr**0.6D0))
              ZoLu(m)    = vonKar*g*blk_ZW*
     &             (Tstar*(1.0D0+cpvir*Q)+cpvir*TairK*Qstar)/
     &             (TairK*Wstar(m)*Wstar(m)*(1.0D0+cpvir*Q)+eps)
              psi_u(m)     = bulk_psiu_coare ( ZoLu(m) )
              logus10 = LOG(blk_ZW*iZoW)
              Wstar(m)   = delW(m)*vonKar/(logus10-psi_u(m))
            ENDDO
            ZoLt      = ZoLu(mb)*blk_ZToZW
            psi_t     = bulk_psit_coare ( ZoLt )
            logts10   = LOG(blk_ZT*iZoT)
            cff       = vonKar/(logts10-psi_t)
            Tstar     = delT*cff
            Qstar     = delQ*cff
              Bf=-g/TairK*Wstar(1)*(Tstar+cpvir*TairK*Qstar)
              if (Bf.gt.0.0D0) then
                cff=blk_beta*(Bf*blk_Zabl)**r3
              else
                cff=0.2D0
              endif
              delW(1)  = SQRT(wspd0*wspd0+cff*cff)
            ddelW = delW
          ENDDO
          hfsen = - blk_Cpa*rhoAir*Wstar(mb)*Tstar
          Hlv   = (2.5008D0 - 0.0023719D0*TseaC)*1.0D+6
          hflat = - Hlv*rhoAir*Wstar(mb)*Qstar
          hflw   = radlw(i,j)
     &           - emiss_lw*rho0i*cpi*SigmaSB*TseaK*TseaK*TseaK*TseaK
          upvel=-1.61D0*Wstar(mb)*Qstar-(1.0D0+1.61D0*Q)*Wstar(mb)
     &                                           *     Tstar/TairK
          hflat=hflat+rhoAir*Hlv*upvel*Q
          hflat=-hflat*rho0i*cpi
          hfsen=-hfsen*rho0i*cpi
          stflx(i,j,itemp)=srflx(i,j)+hflw+hflat+hfsen
          evap=-cp*hflat/Hlv
          stflx(i,j,isalt)=(evap-prate(i,j))*t(i,j,N,nrhs,isalt)
          stflx(i,j,itemp)=stflx(i,j,itemp)*rmask(i,j)
          stflx(i,j,isalt)=stflx(i,j,isalt)*rmask(i,j)
          aer(i,j)  = rhoAir*delW(1)
          Cd        = (Wstar(1)/ddelW(1))**2
          cer(i,j)  = Cd
          shflx_rsw(i,j)=srflx(i,j)
          shflx_lat(i,j)=hflat
          shflx_sen(i,j)=hfsen
          shflx_rlw(i,j)=hflw
        ENDDO
      ENDDO
      do j=jmin,jmax
        do i=imin+1,imax
          cff =0.5D0*(cer(i-1,j)+cer(i,j))
          cff2=0.5D0*(aer(i-1,j)+aer(i,j))
          sustr(i,j)= ( cff*cff2*uwnd(i,j)
     &                             )*rho0i
          sustr(i,j)=sustr(i,j)*umask(i,j)
        enddo
      enddo
      do j=jmin+1,jmax
        do i=imin,imax
          cff =0.5D0*(cer(i,j-1)+cer(i,j))
          cff2=0.5D0*(aer(i,j-1)+aer(i,j))
          svstr(i,j)=( cff*cff2*vwnd(i,j)
     &                                        )*rho0i
          svstr(i,j)=svstr(i,j)*vmask(i,j)
        enddo
      enddo
      return
       
      

      end subroutine Sub_Loop_bulk_flux_tile

      FUNCTION spec_hum (RH,psfc,TairC)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      REAL :: RH
      REAL :: psfc
      REAL :: TairC
      REAL :: spec_hum
                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                           
      REAL :: cff
                                 
      cff=(1.0007D0+3.46D-6*0.01D0*psfc)*6.1121D0*
     &        exp(17.502D0*TairC/(240.97D0+TairC))
      IF (RH.lt.2.0D0) then
         cff=cff*RH
         spec_hum=MvoMa*(cff/(psfc*0.01D0-0.378D0*cff))
      ELSE
         spec_hum=0.001D0*RH
      ENDIF
      END FUNCTION spec_hum
      SUBROUTINE exner_patm_from_tairabs (iexn,pair,q,tairabs,z,psfc)
      use Agrif_Types, only : Agrif_tabvars

      implicit none

      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      REAL, parameter :: g = 9.80665D0
      INTEGER(4), parameter :: Niter = 3
      REAL, intent(out) :: iexn
      REAL, intent(out) :: pair
      REAL, intent(in) :: q
      REAL, intent(in) :: tairabs
      REAL, intent(in) :: z
      REAL, intent(in) :: psfc

                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                        
                                                  
                                            
      REAL :: xm
      REAL :: q_sat
      REAL :: qsat
                                                
                                     
      INTEGER(4) :: iter
                                          
      pair = psfc
      DO Iter = 1, Niter
        q_sat = qsat(tairabs, pair, 1.D0)
        xm    =  mm_dryair + (q/q_sat) * ( mm_water - mm_dryair )
        pair  = psfc * EXP( -g * xm * z / ( r_gas * tairabs ) )
      ENDDO
      iexn =  (pair*ip00)**(-rdocpd)
      return
      END SUBROUTINE exner_patm_from_tairabs
      FUNCTION qsat (TairK, patm, coeff)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      REAL, parameter :: alpw = 60.2227554D0
      REAL, parameter :: betaw = 6822.40088D0
      REAL, parameter :: gamw = 5.13926744D0
      REAL, parameter :: alpi = 32.62117980819471D0
      REAL, parameter :: betai = 6295.421338904806D0
      REAL, parameter :: gami = 0.5631331575423155D0
      REAL :: TairK
      REAL :: patm
      REAL :: coeff
      REAL :: qsat
                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                    
                                                  
                                    
      REAL :: psat
                                                      
                                                      
                                                      
                                                             
                                                             
                                                              
      IF (TairK .LE. CtoK) then
        psat = EXP( alpi - betai/TairK - gami*LOG(TairK) )
      ELSE
        psat = EXP( alpw - betaw/TairK - gamw*LOG(TairK) )
      ENDIF
      psat = coeff * psat
      qsat = (MvoMa*psat)/(patm+(MvoMa-1.0D0)*psat)
      return
      END FUNCTION qsat
      FUNCTION air_visc(TairC)

      use Agrif_Types, only : Agrif_tabvars


      REAL, parameter :: c0 = 1.326D-5
      REAL, parameter :: c1 = 6.542D-3
      REAL, parameter :: c2 = 8.301D-6
      REAL, parameter :: c3 = 4.84D-9
      REAL :: air_visc
                                           
      REAL :: cff
                                            
                                            
                                            
                                           
      cff      = TairC*TairC
      air_visc = c0*(1.D0+c1*TairC+c2*cff-c3*cff*TairC)
      return
      END FUNCTION air_visc
      FUNCTION bulk_psiu_coare (ZoL)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      REAL :: ZoL
      REAL :: bulk_psiu_coare
                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                               
                                   
                                          
      REAL :: chik
      REAL :: psik
                                          
      REAL :: chic
      REAL :: psic
      IF (ZoL <= 0.0D0) then
        chik = (1.0D0-15.0D0*ZoL)**0.25D0
        psik = 2.0D0*LOG(0.5D0*(1.0D0+chik))+LOG(0.5D0*(1.0D0+chik**2))
     &                                  -2.0D0*ATAN(chik)+pis2
        chic = (1.0D0-10.15D0*ZoL)**r3
        psic  = 1.5D0*LOG(r3*(chic**2+chic+1.0D0))
     &          - sqr3*ATAN((2.0D0*chic+1.0D0)/sqr3)+2.D0*pis2osqr3
        bulk_psiu_coare=psic+(psik-psic)/(1.0D0+ZoL**2)
      ELSE
        chic=-MIN(50.0D0,0.35D0*ZoL)
        bulk_psiu_coare=-((1.0D0+ZoL)+0.6667D0*(ZoL-14.28D0)*EXP(chic)+
     &                                                          8.525D0)
      ENDIF
      return
      END FUNCTION bulk_psiu_coare
      FUNCTION bulk_psit_coare (ZoL)

      use Agrif_Types, only : Agrif_tabvars

      implicit none

      REAL, parameter :: CtoK = 273.16D0
      REAL, parameter :: blk_Rgas = 287.0596736665907D0
      REAL, parameter :: blk_Rvap = 461.5249933083879D0
      REAL, parameter :: blk_Cpa = 1004.708857833067D0
      REAL, parameter :: ip00 = 1.D-5
      REAL, parameter :: p00 = 1.D+5
      REAL, parameter :: rdocpd = blk_Rgas/blk_Cpa
      REAL, parameter :: cpdord = blk_Cpa/blk_Rgas
      REAL, parameter :: r_gas = 8.314510D0
      REAL, parameter :: mm_dryair = 28.9644D-3
      REAL, parameter :: mm_water = 18.0153D-3
      REAL, parameter :: cpvir = blk_Rvap/blk_Rgas-1.D0
      REAL, parameter :: MvoMa = mm_water/mm_dryair
      REAL, parameter :: eps = 1.D-8
      REAL, parameter :: r3 = 1.D0/3.D0
      REAL, parameter :: pis2 = 2.D0*ATAN(1.D0)
      REAL, parameter :: sqr3 = SQRT(3.D0)
      REAL, parameter :: pis2osqr3 = pis2/sqr3
      REAL, parameter :: LMOmin = -200.D0
      REAL, parameter :: LMOmax = 0.25D0
      REAL, parameter :: blk_beta = 1.2D0
      REAL, parameter :: blk_Zabl = 600.D0
      REAL, parameter :: dWstar0 = 1.D-6
      REAL, parameter :: dTstar0 = 1.D-6
      REAL, parameter :: dQstar0 = 1.D-9
      REAL, parameter :: emiss_lw = 0.985D0
      REAL, parameter :: SigmaSB = 5.6697D-8
      REAL, parameter :: psurf = 100000.0D0
      REAL :: ZoL
      REAL :: bulk_psit_coare
                   

                                                
                                                           
                                                           
                                                          
                                               
                                               
                                                       
                                                       
                                                 
                                                 
                                                 
                                                               
                                                         
                                            
                                                
                                                      
                                                 
                                                
                                              
                                             
                                            
                                             
                                            
                                            
                                            
                                               
                                                 
                               
                                  

                                               
                                   
                                          
      REAL :: chik
      REAL :: psik
                                          
      REAL :: chic
      REAL :: psic
      IF (ZoL < 0.0D0) THEN
        chik  = (1.0D0-15.0D0*ZoL)**0.25D0
        psik  = 2.0D0*LOG(0.5D0*(1.0D0+chik**2))
        chic  = (1.0D0-34.15D0*ZoL)**r3
        psic  = 1.5D0*LOG((chic**2+chic+1.0D0)*r3)
     &               -sqr3*ATAN((2.0D0*chic+1.0D0)/sqr3)
     &               +2.D0*pis2osqr3
        bulk_psit_coare = psic+(psik-psic)/(1.0D0+ZoL**2)
      ELSE
        chic=-MIN(50.0D0,0.35D0*ZoL)
        bulk_psit_coare = -((1.0D0+2.0D0*ZoL/3.0D0)**1.5D0+
     &            0.6667D0*(ZoL-14.28D0)*EXP(chic)+8.525D0)
      ENDIF
      return
      END FUNCTION bulk_psit_coare
