Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:43 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b2ea1d2979 Jean*0001 module mixed_layer_mod
                0002 
                0003 !
                0004 ! Implementation of mixed layer boundary condition
                0005 !
                0006 
                0007 !use                  fms_mod, only: set_domain, write_version_number, &
                0008 !                                    mpp_pe, mpp_root_pe, error_mesg, FATAL, WARNING
                0009 
                0010 !use                  fms_mod, only: stdlog, check_nml_error, close_file,&
                0011 !                                    open_namelist_file, stdout, file_exist, &
                0012 !                                    read_data, write_data, open_file
                0013 
                0014 use           gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0015 use            constants_mod, only: HLV, PI, RHO_CP, CP_AIR
                0016 
                0017 !use         diag_manager_mod, only: register_diag_field, send_data
                0018 
                0019 !use       time_manager_mod,   only: time_type
                0020 
                0021 !use           transforms_mod, only: get_deg_lat
                0022 
                0023 !use            vert_diff_mod, only: surf_diff_type
                0024 
                0025 implicit none
                0026 private
                0027 !===================================================================================================
                0028 
15388b052e Jean*0029 character(len=128) :: version = '$Id: mixed_layer_mod.F90,v 1.3 2017/08/11 20:48:51 jmc Exp $'
b2ea1d2979 Jean*0030 character(len=128) :: tagname = '$Name:  $'
                0031 character(len=128), parameter :: mod_name='mixed_layer'
                0032 
                0033 !===================================================================================================
                0034 
                0035 public :: mixed_layer_init, mixed_layer, mixed_layer_end
                0036 
                0037 !===================================================================================================
                0038 
                0039 logical :: evaporation = .true.
                0040 real    :: qflux_amp = 0.0
                0041 real    :: qflux_width = 16.0  ! width of qflux region in degrees
                0042 real    :: depth = 40.0
                0043 
                0044 namelist/mixed_layer_nml/ evaporation, qflux_amp, depth, qflux_width
                0045 
                0046 !===================================================================================================
                0047 
                0048 logical :: module_is_initialized =.false.
                0049 logical :: used
                0050 
                0051 integer :: iter
                0052 integer, dimension(4) :: axes
                0053 !integer ::                                                                    &
                0054 !     id_t_surf,            &   ! surface temperature
                0055 !     id_flux_lhe,          &   ! latent heat flux at surface
                0056 !     id_flux_oceanq,       &   ! oceanic Q flux
                0057 !     id_flux_t                 ! sensible heat flux at surface
                0058 
                0059 !real, allocatable, dimension(:,:)   ::                                        &
                0060 !     ocean_qflux,           &   ! Q-flux
                0061 !     rad_lat_2d                 ! latitude in radians
                0062 
                0063 !real, allocatable, dimension(:)   :: deg_lat
                0064 
                0065 !real, allocatable, dimension(:,:)   ::                                        &
                0066 !     gamma_t,               &   ! Used to calculate the implicit
                0067 !     gamma_q,               &   ! correction to the diffusion in
                0068 !     fn_t,                  &   ! the lowest layer
                0069 !     fn_q,                  &   !
                0070 !     en_t,                  &   !
                0071 !     en_q,                  &   !
                0072 !     alpha_t,               &   !
                0073 !     alpha_q,               &   !
                0074 !     alpha_lw,              &   !
                0075 !     beta_t,                &   !
                0076 !     beta_q,                &   !
                0077 !     beta_lw,               &   !
                0078 !     t_surf_dependence,     &   !
                0079 !     corrected_flux,        &   !
                0080 !     eff_heat_capacity,     &   ! Effective heat capacity
                0081 !     delta_t_surf               ! Increment in surface temperature
                0082 
                0083 real inv_cp_air
                0084 
                0085 !===================================================================================================
                0086 contains
                0087 !===================================================================================================
                0088 
                0089 !subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, axes,  &
                0090 !                            rad_lat_2d, ocean_qflux,                   &
c9694dc201 Jean*0091 subroutine mixed_layer_init( is, ie, js, je, num_levels,         axes,  &
                0092                              Time, cst_mxlDepth, myThid )
b2ea1d2979 Jean*0093 
c9694dc201 Jean*0094 integer, intent(in) :: is, ie, js, je, num_levels
b2ea1d2979 Jean*0095 !real, intent(inout), dimension(:,:) :: t_surf
                0096 integer, intent(in), dimension(4) :: axes
c9694dc201 Jean*0097 !type(time_type), intent(in)       :: Time
                0098 real, intent(in)                   :: Time
                0099 real, intent(out)                  :: cst_mxlDepth
b2ea1d2979 Jean*0100 integer, intent(in)                :: myThid
                0101 !real, intent(in), dimension(:,:)  :: rad_lat_2d
                0102 !real, intent(out), dimension(:,:) :: ocean_qflux
                0103 
c9694dc201 Jean*0104 !integer :: j
                0105 !real    :: rad_qwidth
                0106 !integer:: ierr, io, unit
b2ea1d2979 Jean*0107 integer         :: iUnit
                0108 CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0109 
                0110 if(module_is_initialized) return
                0111 
                0112 ! read namelist and copy to logfile
                0113 !    _BARRIER
                0114 !    _BEGIN_MASTER(myThid)
                0115      CALL BARRIER(myThid)
                0116      IF ( myThid.EQ.1 ) THEN
                0117 
                0118      WRITE(msgBuf,'(A)') 'MIXED_LAYER_INIT: opening data.atm_gray'
                0119      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0120      CALL OPEN_COPY_DATA_FILE(                                      &
c9694dc201 Jean*0121                            'data.atm_gray', 'MIXED_LAYER_INIT',     &
b2ea1d2979 Jean*0122                            iUnit,                                   &
                0123                            myThid )
                0124 !    Read parameters from open data file
                0125      READ(UNIT=iUnit,NML=mixed_layer_nml)
                0126      WRITE(msgBuf,'(A)')                                            &
                0127           'MIXED_LAYER_INIT: finished reading data.atm_gray'
                0128      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0129 !    Close the open data file
15388b052e Jean*0130 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0131      CLOSE(iUnit)
15388b052e Jean*0132 #else
                0133      CLOSE(iUnit,STATUS='DELETE')
                0134 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0135 
                0136 !call write_version_number(version, tagname)
                0137 
                0138 !unit = open_namelist_file ()
                0139 !ierr=1
                0140 !do while (ierr /= 0)
                0141 !  read  (unit, nml=mixed_layer_nml, iostat=io, end=10)
                0142 !  ierr = check_nml_error (io, 'mixed_layer_nml')
                0143 !enddo
                0144 !10 call close_file (unit)
                0145 
                0146 !if ( mpp_pe() == mpp_root_pe() )   write (stdlog(), nml=mixed_layer_nml)
                0147 
                0148 !allocate(rad_lat_2d              (is:ie, js:je))
                0149 !allocate(ocean_qflux             (is:ie, js:je))
                0150 
                0151 !allocate(deg_lat                 (js:je))
                0152 
                0153 !allocate(gamma_t                 (is:ie, js:je))
                0154 !allocate(gamma_q                 (is:ie, js:je))
                0155 !allocate(en_t                    (is:ie, js:je))
                0156 !allocate(en_q                    (is:ie, js:je))
                0157 !allocate(fn_t                    (is:ie, js:je))
                0158 !allocate(fn_q                    (is:ie, js:je))
                0159 !allocate(alpha_t                 (is:ie, js:je))
                0160 !allocate(alpha_q                 (is:ie, js:je))
                0161 !allocate(alpha_lw                (is:ie, js:je))
                0162 !allocate(beta_t                  (is:ie, js:je))
                0163 !allocate(beta_q                  (is:ie, js:je))
                0164 !allocate(beta_lw                 (is:ie, js:je))
                0165 !allocate(delta_t_surf            (is:ie, js:je))
                0166 !allocate(eff_heat_capacity       (is:ie, js:je))
                0167 !allocate(corrected_flux          (is:ie, js:je))
                0168 !allocate(t_surf_dependence       (is:ie, js:je))
                0169 
                0170 inv_cp_air = 1.0 / CP_AIR
                0171 
                0172 module_is_initialized = .true.
                0173 
                0174      ENDIF
                0175      CALL BARRIER(myThid)
c9694dc201 Jean*0176      cst_mxlDepth = depth
b2ea1d2979 Jean*0177 
                0178 !see if restart file exists for the surface temperature
                0179 !
                0180 !if (file_exist('INPUT/mixed_layer.res')) then
                0181 !         unit = open_file (file='INPUT/mixed_layer.res', &
                0182 !                           form='native', action='read')
                0183 !         call read_data (unit, t_surf)
                0184 !         call close_file (unit)
                0185 !else if (file_exist('INPUT/swamp.res')) then
                0186 !         unit = open_file (file='INPUT/swamp.res', &
                0187 !                           form='native', action='read')
                0188 !         call read_data (unit, t_surf)
                0189 !         call close_file (unit)
                0190 !  call error_mesg('mixed_layer','mixed_layer restart file not found, using swamp restart file', WARNING)
                0191 !else
                0192 !  call error_mesg('mixed_layer','mixed_layer restart file not found', WARNING)
                0193 !endif
                0194 
                0195 !id_t_surf = register_diag_field(mod_name, 't_surf',        &
                0196 !                                axes(1:2), Time, 'surface temperature','K')
                0197 !id_flux_t = register_diag_field(mod_name, 'flux_t',        &
                0198 !                                axes(1:2), Time, 'sensible heat flux up at surface','watts/m2')
                0199 !id_flux_lhe = register_diag_field(mod_name, 'flux_lhe',        &
                0200 !                                 axes(1:2), Time, 'latent heat flux up at surface','watts/m2')
                0201 !id_flux_oceanq = register_diag_field(mod_name, 'flux_oceanq',        &
                0202 !                                 axes(1:2), Time, 'oceanic Q-flux','watts/m2')
                0203 
                0204 ! latitude will be needed for oceanic q flux
                0205 !call get_deg_lat(deg_lat)
                0206 !do j=js,je
                0207 !  rad_lat_2d(:,j) = deg_lat(j)*PI/180.
                0208 !enddo
                0209 
                0210 ! calculate ocean Q flux
                0211 !rad_qwidth = qflux_width*PI/180.
                0212 !ocean_qflux = qflux_amp*(1-2.*rad_lat_2d**2/rad_qwidth**2) * &
                0213 !        exp(- ((rad_lat_2d)**2/(rad_qwidth)**2))
                0214 
                0215 return
                0216 end subroutine mixed_layer_init
                0217 
                0218 !===================================================================================================
                0219 
                0220 subroutine mixed_layer (                                &
                0221      Time,                                              &
                0222      t_surf,                                            &
                0223      flux_t,                                            &
                0224      flux_q,                                            &
                0225      flux_r,                                            &
                0226      dt,                                                &
                0227      net_surf_sw_down,                                  &
                0228      surf_lw_down,                                      &
                0229 !    Tri_surf,                                          &
                0230                 tri_surf_dtmass,                               &
                0231                 tri_surf_dflux_t, tri_surf_dflux_q,            &
                0232                 tri_surf_delta_t, tri_surf_delta_q,            &
                0233      dhdt_surf,                                         &
                0234      dedt_surf,                                         &
                0235      dedq_surf,                                         &
                0236      drdt_surf,                                         &
                0237      dhdt_atm,                                          &
                0238      dedq_atm,                                          &
c9694dc201 Jean*0239      ocean_qflux, mixLayDepth,                          &
b2ea1d2979 Jean*0240      delta_t_surf, &               ! Increment in surface temperature
                0241      myThid )
                0242 
                0243 ! ---- arguments -----------------------------------------------------------
                0244 !type(time_type), intent(in)       :: Time
                0245 real, intent(in)                  :: Time
                0246 real, intent(in),  dimension(:,:) :: &
                0247      net_surf_sw_down, surf_lw_down
                0248 real, intent(inout), dimension(:,:) :: &
                0249      flux_t,    flux_q,     flux_r
                0250 real, intent(inout), dimension(:,:) :: t_surf
                0251 real, intent(in), dimension(:,:) :: &
                0252    dhdt_surf, dedt_surf, dedq_surf, &
                0253    drdt_surf, dhdt_atm, dedq_atm
                0254 real, intent(in) :: dt
                0255 !type(surf_diff_type), intent(inout) :: Tri_surf
                0256 real, intent(inout), dimension(:,:)  :: tri_surf_dtmass
                0257 real, intent(inout), dimension(:,:)  :: tri_surf_dflux_t, tri_surf_dflux_q
                0258 real, intent(inout), dimension(:,:)  :: tri_surf_delta_t, tri_surf_delta_q
                0259 real, intent(in),    dimension(:,:)  :: ocean_qflux
c9694dc201 Jean*0260 real, intent(in),    dimension(:,:)  :: mixLayDepth
b2ea1d2979 Jean*0261 real, intent(out),   dimension(:,:)  :: delta_t_surf
                0262 integer, intent(in) :: myThid
                0263 
                0264 ! ---- local variables --------------------------------------------------
                0265 real, allocatable, dimension(:,:)   ::                                        &
                0266      gamma_t,               &   ! Used to calculate the implicit
                0267      gamma_q,               &   ! correction to the diffusion in
                0268      fn_t,                  &   ! the lowest layer
                0269      fn_q,                  &   !
                0270      en_t,                  &   !
                0271      en_q,                  &   !
                0272      alpha_t,               &   !
                0273      alpha_q,               &   !
                0274      alpha_lw,              &   !
                0275      beta_t,                &   !
                0276      beta_q,                &   !
                0277      beta_lw,               &   !
                0278      t_surf_dependence,     &   !
                0279      corrected_flux,        &   !
                0280      eff_heat_capacity          ! Effective heat capacity
                0281 !    eff_heat_capacity,     &   ! Effective heat capacity
                0282 !    delta_t_surf               ! Increment in surface temperature
                0283 integer :: im, jm
                0284 !----------------------------------------------------------------
                0285 
                0286 im = size(t_surf,1)
                0287 jm = size(t_surf,2)
                0288 allocate(gamma_t            (im,jm) )
                0289 allocate(gamma_q            (im,jm) )
                0290 allocate(en_t               (im,jm) )
                0291 allocate(en_q               (im,jm) )
                0292 allocate(fn_t               (im,jm) )
                0293 allocate(fn_q               (im,jm) )
                0294 allocate(alpha_t            (im,jm) )
                0295 allocate(alpha_q            (im,jm) )
                0296 allocate(alpha_lw           (im,jm) )
                0297 allocate(beta_t             (im,jm) )
                0298 allocate(beta_q             (im,jm) )
                0299 allocate(beta_lw            (im,jm) )
                0300 allocate(eff_heat_capacity  (im,jm) )
                0301 allocate(corrected_flux     (im,jm) )
                0302 allocate(t_surf_dependence  (im,jm) )
                0303 !allocate(delta_t_surf       (im,jm) )
                0304 
                0305 if(.not.module_is_initialized) then
                0306 ! call error_mesg('mixed_layer','mixed_layer module is not initialized',FATAL)
                0307   call PRINT_ERROR('mixed_layer: mixed_layer module is not initialized',myThid)
                0308   STOP 'ABNORMAL END: S/R MIXED_LAYER'
                0309 endif
                0310 
                0311 ! Need to calculate the implicit changes to the lowest level delta_q and delta_t
                0312 ! - see the discussion in vert_diff.tech.ps
                0313 
                0314 ! Care is needed to differentiate between the sensible heat flux and the
                0315 ! diffusive flux of temperature
                0316 
                0317 gamma_t = 1.0 / (1.0 - Tri_surf_dtmass * (Tri_surf_dflux_t + dhdt_atm * inv_cp_air))
                0318 gamma_q = 1.0 / (1.0 - Tri_surf_dtmass * (Tri_surf_dflux_q + dedq_atm))
                0319 
                0320 fn_t = gamma_t * (Tri_surf_delta_t + Tri_surf_dtmass * flux_t * inv_cp_air)
                0321 fn_q = gamma_q * (Tri_surf_delta_q + Tri_surf_dtmass * flux_q)
                0322 
                0323 en_t = gamma_t * Tri_surf_dtmass * dhdt_surf * inv_cp_air
                0324 en_q = gamma_q * Tri_surf_dtmass * dedt_surf
                0325 
                0326 !
                0327 ! Note flux_sw doesn't depend on surface or lowest layer values
                0328 ! Note drdt_atm is not used - should be fixed
                0329 !
                0330 alpha_t = flux_t * inv_cp_air + dhdt_atm * inv_cp_air * fn_t
                0331 alpha_q = flux_q + dedq_atm * fn_q
                0332 alpha_lw = flux_r
                0333 
                0334 beta_t = dhdt_surf * inv_cp_air + dhdt_atm * inv_cp_air * en_t
                0335 beta_q = dedt_surf + dedq_atm * en_q
                0336 beta_lw = drdt_surf
                0337 
                0338 !
                0339 ! Implement mixed layer surface boundary condition
                0340 !
                0341 corrected_flux = - net_surf_sw_down - surf_lw_down + alpha_t * CP_AIR + alpha_lw + ocean_qflux
                0342 t_surf_dependence = beta_t * CP_AIR + beta_lw
                0343 
                0344 if (evaporation) then
                0345   corrected_flux = corrected_flux + alpha_q * HLV
                0346   t_surf_dependence = t_surf_dependence + beta_q * HLV
                0347 endif
                0348 
                0349 !
                0350 ! Now update the mixed layer surface temperature using an implicit step
                0351 !
c9694dc201 Jean*0352 !eff_heat_capacity = depth * RHO_CP + t_surf_dependence * dt
                0353 eff_heat_capacity = mixLayDepth * RHO_CP + t_surf_dependence * dt
b2ea1d2979 Jean*0354 
                0355 if (any(eff_heat_capacity .eq. 0.0))  then
                0356   write(*,*) 'mixed_layer: error', eff_heat_capacity
                0357 ! call error_mesg('mixed_layer', 'Avoiding division by zero',fatal)
                0358   call PRINT_ERROR('mixed_layer: Avoiding division by zero',myThid)
                0359   STOP 'ABNORMAL END: S/R MIXED_LAYER'
                0360 end if
                0361 
                0362 delta_t_surf = - corrected_flux  * dt / eff_heat_capacity
                0363 
                0364 t_surf = t_surf + delta_t_surf
                0365 
                0366 ! Finally calculate the increments for the lowest atmospheric layer
                0367 !
                0368 Tri_surf_delta_t = fn_t + en_t * delta_t_surf
                0369 Tri_surf_delta_q = fn_q + en_q * delta_t_surf
                0370 
                0371 ! Note:
                0372 ! When using an implicit step there is not a clearly defined flux for a given timestep
                0373 
                0374 ! Final flux values (match mixed layer heat content variations)
                0375   flux_t = ( alpha_t  + delta_t_surf*beta_t )*CP_AIR
                0376   flux_r =   alpha_lw + delta_t_surf*beta_lw
                0377   flux_q =   alpha_q  + delta_t_surf*beta_q
                0378 
                0379 !if(id_t_surf > 0) used = send_data(id_t_surf, t_surf, Time)
                0380 !if(id_flux_t > 0) used = send_data(id_flux_t, flux_t, Time)
                0381 !if(id_flux_lhe > 0) used = send_data(id_flux_lhe, HLV * flux_q, Time)
                0382 !if(id_flux_oceanq > 0)   used = send_data(id_flux_oceanq, ocean_qflux, Time)
                0383 
                0384 deallocate(gamma_t            )
                0385 deallocate(gamma_q            )
                0386 deallocate(en_t               )
                0387 deallocate(en_q               )
                0388 deallocate(fn_t               )
                0389 deallocate(fn_q               )
                0390 deallocate(alpha_t            )
                0391 deallocate(alpha_q            )
                0392 deallocate(alpha_lw           )
                0393 deallocate(beta_t             )
                0394 deallocate(beta_q             )
                0395 deallocate(beta_lw            )
                0396 deallocate(eff_heat_capacity  )
                0397 deallocate(corrected_flux     )
                0398 deallocate(t_surf_dependence  )
                0399 !deallocate(delta_t_surf       )
                0400 
                0401 end subroutine mixed_layer
                0402 
                0403 !===================================================================================================
                0404 
                0405 subroutine mixed_layer_end(t_surf,myThid)
                0406 
                0407 real, intent(inout), dimension(:,:) :: t_surf
                0408 integer, intent(in)                :: myThid
                0409 
c9694dc201 Jean*0410 !integer:: unit
b2ea1d2979 Jean*0411 
                0412 if(.not.module_is_initialized) return
                0413 
                0414 ! write a restart file for the surface temperature
                0415 !unit = open_file ('RESTART/mixed_layer.res', &
                0416 !                  form='native', action='write')
                0417 !call write_data (unit, t_surf)
                0418 !call close_File (unit)
                0419 
                0420 module_is_initialized = .false.
                0421 
                0422 end subroutine mixed_layer_end
                0423 
                0424 !===================================================================================================
                0425 
                0426 end module mixed_layer_mod
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!