<commit_before>
module neural_convection_mod


!----------------------------------------------------------------------
use            fms_mod, only:  file_exist, error_mesg, open_file,  &
                               check_nml_error, mpp_pe, FATAL,  &
                               close_file
use simple_sat_vapor_pres_mod, only:  escomp, descomp
use      constants_mod, only:  HLv,HLs,Cp_air,Grav,rdgas,rvgas, &
                               kappa
use netcdf

implicit none
private
!---------------------------------------------------------------------
!  ---- public interfaces ----

   public  neural_convection, neural_convection_init, check

!-----------------------------------------------------------------------
!   ---- version number ----

 character(len=128) :: version = '$Id: neural_convection.f90,v 1 2016/05/24 00:00:00 fms Exp $'
 character(len=128) :: tag = '$Name: fez $'

!-----------------------------------------------------------------------
!   ---- local/private data ----

    logical :: do_init=.true.

!-----------------------------------------------------------------------
!   ---- namelist ----

character(len=128) :: neural_filename = 'INPUT/neural_weights_v3.nc'

namelist /neural_convection_nml/ neural_filename

contains

!#######################################################################

   subroutine neural_convection (dt, tin, qin, pfull, phalf,  &
                                  rain, tdel, qdel,  &
                                   r_w1, r_w2, r_b1, r_b2, &
                               xscale_mean,xscale_stnd,yscale_absmax)

!-----------------------------------------------------------------------
!
!                     Betts-Miller Convection Scheme
!
!-----------------------------------------------------------------------
!
!   input:  dt       time step in seconds
!           tin      temperature at full model levels
!           qin      specific humidity of water vapor at full
!                      model levels
!           pfull    pressure at full model levels
!           phalf    pressure at half (interface) model levels
!
!  output:  rain     total precipitation (kg/m2)
!           tdel     temperature tendency at full model levels
!           qdel     specific humidity tendency (of water vapor) at
!                      full model levels
!
!-----------------------------------------------------------------------
!--------------------- interface arguments -----------------------------

   real   , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf
   real   , intent(in)                    :: dt
   real   , intent(out), dimension(:,:)   :: rain
   real   , intent(out), dimension(:,:,:) :: tdel, qdel
!-----------------------------------------------------------------------
!---------------------- local data -------------------------------------

   real,dimension(size(tin,1),size(tin,2))             :: precip
   real,dimension(size(tin,3)/2)                         :: qpc, tpc
   real,dimension(size(tin,3))                       :: features, targets
   real                                                :: deltak, qnew_tmp
 integer  i, j, k, ix, jx, kx, ktop, kx2

      real, intent(in), dimension(size(tin,3), 100)   :: r_w1 
      real, intent(in), dimension(100)                :: r_b1 
      real, intent(in), dimension(100, size(tin,3))   :: r_w2
      real, intent(in), dimension(size(tin,3))        :: r_b2
      real, intent(in), dimension(size(tin,3))        :: xscale_mean
      real, intent(in), dimension(size(tin,3))        :: xscale_stnd
      real, intent(in), dimension(size(tin,3))        :: yscale_absmax

      real, dimension(100) :: z1, z2
!     real, dimension(300) :: a1, a2
!      real, dimension(2)   :: softmax
!      real                 :: prob
!-----------------------------------------------------------------------
!     computation of precipitation by betts-miller scheme
!-----------------------------------------------------------------------

      if (do_init) call error_mesg ('neural_convection',  &
                         'neural_convection_init has not been called.', FATAL)
      ix=size(tin,1)
      jx=size(tin,2)
      kx=size(tin,3)
      kx2 = kx / 2
       do i=1,ix
          do j=1,jx
! Initialize variables
             precip(i,j) = 0.
             features = 0.
!             a1 = 0.
!             a2 = 0.
             qpc = 0.
             tpc = 0.
             z1 = 0.
             z2 = 0.
!             softmax = 0.
!             prob = 0.
             targets = 0.
! Temperature and humidity profiles
             tpc = tin(i,j,kx2:kx)
             qpc = qin(i,j,kx2:kx)

! Combine tpc and rpc into a vector of levels (2*N_lev)
             features(1:kx2) = tpc
             features(kx2+1:kx) = qpc
! Scale inputs
! See http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html#sklearn.preprocessing.MinMaxScaler.get_params for more info about scaling
!             features = (features - xscale_min) / (xscale_max - xscale_min)
!             features = features * 2.0 - 1.0 
             features = (features - xscale_mean) / xscale_stnd

! Classify data
!             a1 = matmul(features,  c_w1) + c_b1 
!             where (a1 .lt. 0.0)  a1 = 0.0
!             a2 = matmul(a1      ,  c_w2) + c_b2
!             where (a2 .lt. 0.0)  a2 = 0.0
!             softmax = matmul(a2,   c_w3) + c_b3
! Use exponentials to calculate the probabilities via the softmax function
!             softmax = exp(softmax)
!             prob = softmax(1) / ( softmax(1) + softmax(2) )
! If classifier says convection is occurring, run model with regressor
!            if (prob .gt. 0.5) then
!                 z1 = matmul(features,  r_w1) + r_b1
!                 where (z1 .lt. 0.0)  z1 = 0.0
!                 z2 = matmul(z1,        r_w2) + r_b2
!                 where (z2 .lt. 0.0)  z2 = 0.0
!                 targets = matmul(z2, r_w3) + r_b3
!             else
!                 targets = 0.
!             endif

! Apply trained regressor network to data using rectifier activation function
! forward prop to hidden layer
             z1 = matmul(features,  r_w1) + r_b1
! rectifier
             where (z1 .lt. 0.0)  z1 = 0.0
! forward prop to output layer
             targets = matmul(z1, r_w2) + r_b2
! Inverse scale outputs
             targets = targets * yscale_absmax
! Separate out targets into heating and moistening tendencies
             tdel(i,j,kx2:kx) = targets(1:kx2)
             qdel(i,j,kx2:kx) = targets(kx2 + 1:kx)
! Correct units (nn trained on K/day and g/kg/day rather than K and kg/kg)
             tdel(i,j,:) = tdel(i,j,:) / 3600. / 24.
             qdel(i,j,:) = qdel(i,j,:) / 3600. / 24. / 1000.
! If any humidities would become negative set them to zero (and recalc precip)
             do k=1, kx2
                 qnew_tmp = 0.0
                 qnew_tmp = qpc(k) + qdel(i,j,kx2 + k)
                 if ( qnew_tmp .lt. 0.0 ) then
                     ! Dry out level, but don't let it go negative<commit_msg>Accidentally introduced dt bug on prevous commit. Fixed that now. Also updated to deal with 16 levels (rather than 15) and made the code work more generally for levels other than 15<commit_after>
module neural_convection_mod


!----------------------------------------------------------------------
use            fms_mod, only:  file_exist, error_mesg, open_file,  &
                               check_nml_error, mpp_pe, FATAL,  &
                               close_file
use simple_sat_vapor_pres_mod, only:  escomp, descomp
use      constants_mod, only:  HLv,HLs,Cp_air,Grav,rdgas,rvgas, &
                               kappa
use netcdf

implicit none
private
!---------------------------------------------------------------------
!  ---- public interfaces ----

   public  neural_convection, neural_convection_init, check

!-----------------------------------------------------------------------
!   ---- version number ----

 character(len=128) :: version = '$Id: neural_convection.f90,v 1 2016/05/24 00:00:00 fms Exp $'
 character(len=128) :: tag = '$Name: fez $'

!-----------------------------------------------------------------------
!   ---- local/private data ----

    logical :: do_init=.true.

!-----------------------------------------------------------------------
!   ---- namelist ----

character(len=128) :: neural_filename = 'INPUT/neural_weights_v3.nc'

namelist /neural_convection_nml/ neural_filename

contains

!#######################################################################

   subroutine neural_convection (dt, tin, qin, pfull, phalf,  &
                                  rain, tdel, qdel,  &
                                   r_w1, r_w2, r_b1, r_b2, &
                               xscale_mean,xscale_stnd,yscale_absmax)

!-----------------------------------------------------------------------
!
!                     Betts-Miller Convection Scheme
!
!-----------------------------------------------------------------------
!
!   input:  dt       time step in seconds
!           tin      temperature at full model levels
!           qin      specific humidity of water vapor at full
!                      model levels
!           pfull    pressure at full model levels
!           phalf    pressure at half (interface) model levels
!
!  output:  rain     total precipitation (kg/m2)
!           tdel     temperature tendency at full model levels
!           qdel     specific humidity tendency (of water vapor) at
!                      full model levels
!
!-----------------------------------------------------------------------
!--------------------- interface arguments -----------------------------

   real   , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf
   real   , intent(in)                    :: dt
   real   , intent(out), dimension(:,:)   :: rain
   real   , intent(out), dimension(:,:,:) :: tdel, qdel
!-----------------------------------------------------------------------
!---------------------- local data -------------------------------------

   real,dimension(size(tin,1),size(tin,2))             :: precip
   real,dimension(32/2)                         :: qpc, tpc
   real,dimension(32)                       :: features, targets
   real                                                :: deltak, qnew_tmp
 integer  i, j, k, ix, jx, kx, ktop, kx2, kx2ind

      real, intent(in), dimension(32, 100)   :: r_w1 
      real, intent(in), dimension(100)       :: r_b1 
      real, intent(in), dimension(100, 32)   :: r_w2
      real, intent(in), dimension(32)        :: r_b2
      real, intent(in), dimension(32)        :: xscale_mean
      real, intent(in), dimension(32)        :: xscale_stnd
      real, intent(in), dimension(32)        :: yscale_absmax

      real, dimension(100) :: z1, z2
!     real, dimension(300) :: a1, a2
!      real, dimension(2)   :: softmax
!      real                 :: prob
!-----------------------------------------------------------------------
!     computation of precipitation by betts-miller scheme
!-----------------------------------------------------------------------

      if (do_init) call error_mesg ('neural_convection',  &
                         'neural_convection_init has not been called.', FATAL)
      ix=size(tin,1)
      jx=size(tin,2)
      kx=size(tin,3)
      kx2 = size(r_b2)/2  ! Total number of levels the NN uses
      kx2ind = kx - kx2 + 1  ! Index value to grab data for NN from full profile
       do i=1,ix
          do j=1,jx
! Initialize variables
             precip(i,j) = 0.
             features = 0.
!             a1 = 0.
!             a2 = 0.
             qpc = 0.
             tpc = 0.
             z1 = 0.
             z2 = 0.
!             softmax = 0.
!             prob = 0.
             targets = 0.
! Temperature and humidity profiles
             tpc = tin(i,j,kx2ind:kx)
             qpc = qin(i,j,kx2ind:kx)

! Combine tpc and rpc into a vector of levels (2*N_lev)
             features(1:kx2) = tpc
             features(kx2+1:2*kx2) = qpc
! Scale inputs
! See http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html#sklearn.preprocessing.MinMaxScaler.get_params for more info about scaling
!             features = (features - xscale_min) / (xscale_max - xscale_min)
!             features = features * 2.0 - 1.0 
             features = (features - xscale_mean) / xscale_stnd

! Classify data
!             a1 = matmul(features,  c_w1) + c_b1 
!             where (a1 .lt. 0.0)  a1 = 0.0
!             a2 = matmul(a1      ,  c_w2) + c_b2
!             where (a2 .lt. 0.0)  a2 = 0.0
!             softmax = matmul(a2,   c_w3) + c_b3
! Use exponentials to calculate the probabilities via the softmax function
!             softmax = exp(softmax)
!             prob = softmax(1) / ( softmax(1) + softmax(2) )
! If classifier says convection is occurring, run model with regressor
!            if (prob .gt. 0.5) then
!                 z1 = matmul(features,  r_w1) + r_b1
!                 where (z1 .lt. 0.0)  z1 = 0.0
!                 z2 = matmul(z1,        r_w2) + r_b2
!                 where (z2 .lt. 0.0)  z2 = 0.0
!                 targets = matmul(z2, r_w3) + r_b3
!             else
!                 targets = 0.
!             endif

! Apply trained regressor network to data using rectifier activation function
! forward prop to hidden layer
             z1 = matmul(features,  r_w1) + r_b1
! rectifier
             where (z1 .lt. 0.0)  z1 = 0.0
! forward prop to output layer
             targets = matmul(z1, r_w2) + r_b2
! Inverse scale outputs
             targets = targets * yscale_absmax
! Separate out targets into heating and moistening tendencies
             tdel(i,j,kx2ind:kx) = targets(1:kx2)
             qdel(i,j,kx2ind:kx) = targets(kx2 + 1:2*kx2)
! Correct units (nn trained on K/day and g/kg/day rather than K and kg/kg)
             tdel(i,j,:) = tdel(i,j,:) / 3600. / 24.
             qdel(i,j,:) = qdel(i,j,:) / 3600. / 24. / 1000.
! Predicted temperature and humidity tendencies, but actually want the delta T and q
! So multiply each by the time step
             tdel(i,j,:) = tdel(i,j,:) * dt
             qdel(i,j,:) = qdel(i,j,:) * dt
! If any humidities would become negative set them to zero (and recalc precip)
             do k=1, kx2
                 qnew_tmp = 0.0
                 qnew_tmp = qpc(k) + qdel(i,j,kx2 + k)
                 if ( qnew_tmp .lt. 0.0 ) then
                     ! Dry out level, but don't let it go negative
{'fraction_non_alphanumeric': 0.227903125646864, 'fraction_numerical': 0.02808252259711585, 'mean_word_length': 1.2453911696359412, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5799052', 'n_tokens_mistral': 4515, 'n_tokens_neox': 4182, 'n_words': 1107}