48 private :: atmos_refstate_generate_isa
49 private :: atmos_refstate_generate_uniform
50 private :: atmos_refstate_generate_zero
56 character(len=H_LONG),
private :: atmos_refstate_in_basename =
''
57 logical,
private :: atmos_refstate_in_check_coordinates = .true.
58 character(len=H_LONG),
private :: atmos_refstate_out_basename =
''
59 character(len=H_MID),
private :: atmos_refstate_out_title =
'SCALE-RM RefState'
60 character(len=H_SHORT),
private :: atmos_refstate_out_dtype =
'DEFAULT'
62 character(len=H_SHORT),
private :: atmos_refstate_type =
'UNIFORM'
63 real(
rp),
private :: atmos_refstate_temp_sfc = 300.0_rp
64 real(
rp),
private :: atmos_refstate_rh = 0.0_rp
65 real(
rp),
private :: atmos_refstate_pott_uniform = 300.0_rp
66 real(
dp),
private :: atmos_refstate_update_dt = -1.0_dp
68 real(
dp),
private :: last_updated
70 real(
rp),
private,
allocatable :: atmos_refstate1d_pres(:)
71 real(
rp),
private,
allocatable :: atmos_refstate1d_temp(:)
72 real(
rp),
private,
allocatable :: atmos_refstate1d_dens(:)
73 real(
rp),
private,
allocatable :: atmos_refstate1d_pott(:)
74 real(
rp),
private,
allocatable :: atmos_refstate1d_qv (:)
81 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
82 CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
88 integer,
intent(in) :: ka, ks, ke
89 integer,
intent(in) :: ia, is, ie
90 integer,
intent(in) :: ja, js, je
92 real(
rp),
intent(in) :: cz ( ka)
93 real(
rp),
intent(in) :: fz (0:ka)
94 real(
rp),
intent(in) :: real_cz ( ka,ia,ja)
95 real(
rp),
intent(in) :: real_fz (0:ka,ia,ja)
96 real(
rp),
intent(in) :: real_phi( ka,ia,ja)
98 namelist / param_atmos_refstate / &
99 atmos_refstate_in_basename, &
100 atmos_refstate_out_basename, &
101 atmos_refstate_out_title, &
102 atmos_refstate_out_dtype, &
103 atmos_refstate_type, &
104 atmos_refstate_temp_sfc, &
106 atmos_refstate_pott_uniform, &
107 atmos_refstate_update_dt
113 log_info(
"ATMOS_REFSTATE_setup",*)
'Setup'
117 read(
io_fid_conf,nml=param_atmos_refstate,iostat=ierr)
119 log_info(
"ATMOS_REFSTATE_setup",*)
'Not found namelist. Default used.'
120 elseif( ierr > 0 )
then
121 log_error(
"ATMOS_REFSTATE_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_REFSTATE. Check!'
124 log_nml(param_atmos_refstate)
137 allocate( atmos_refstate1d_pres(ka) )
138 allocate( atmos_refstate1d_temp(ka) )
139 allocate( atmos_refstate1d_dens(ka) )
140 allocate( atmos_refstate1d_pott(ka) )
141 allocate( atmos_refstate1d_qv(ka) )
142 atmos_refstate1d_pres(:) = undef
143 atmos_refstate1d_temp(:) = undef
144 atmos_refstate1d_dens(:) = undef
145 atmos_refstate1d_pott(:) = undef
146 atmos_refstate1d_qv(:) = undef
149 log_info(
"ATMOS_REFSTATE_setup",*)
'Reference state settings'
150 if ( atmos_refstate_in_basename /=
'' )
then
151 log_info_cont(*)
'Input file of reference state : ', trim(atmos_refstate_in_basename)
153 log_info_cont(*)
'Input file of reference state : Nothing, generate internally'
156 if ( atmos_refstate_out_basename /=
'' )
then
157 log_info_cont(*)
'Output file of reference state : ', trim(atmos_refstate_out_basename)
159 log_info_cont(*)
'Output file of reference state : No output'
163 if ( atmos_refstate_in_basename /=
'' )
then
167 real_cz(:,:,:), real_fz(:,:,:), &
171 if ( atmos_refstate_type ==
'ISA' )
then
173 log_info_cont(*)
'Reference type : ISA'
174 log_info_cont(*)
'Surface temperature [K] : ', atmos_refstate_temp_sfc
175 log_info_cont(*)
'Surface & environment RH [%] : ', atmos_refstate_rh
177 call atmos_refstate_generate_isa( ka, ks, ke, ia, is, ie, ja, js, je, &
179 real_cz(:,:,:), real_fz(:,:,:), &
182 elseif ( atmos_refstate_type ==
'UNIFORM' )
then
184 log_info_cont(*)
'Reference type : UNIFORM POTT'
185 log_info_cont(*)
'Potential temperature : ', atmos_refstate_pott_uniform
187 call atmos_refstate_generate_uniform( ka, ks, ke, ia, is, ie, ja, js, je, &
189 real_cz(:,:,:), real_fz(:,:,:), &
192 elseif ( atmos_refstate_type ==
'ZERO' )
then
194 log_info_cont(*)
'Reference type : ZERO'
196 call atmos_refstate_generate_zero( ka, ia, ja )
198 elseif ( atmos_refstate_type ==
'INIT' )
then
200 if ( atmos_refstate_update_dt >= 0.0_rp )
then
204 log_info_cont(*)
'Reference type : Generate from initial data'
206 log_info_cont(*)
'Update interval [sec] : ', atmos_refstate_update_dt
209 log_error(
"ATMOS_REFSTATE_setup",*)
'ATMOS_REFSTATE_TYPE must be "ISA", "UNIFORM", "ZERO", or "INIT". Check! : ', trim(atmos_refstate_type)
221 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
222 CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
225 file_cartesc_check_coordinates, &
232 integer,
intent(in) :: ka, ks, ke
233 integer,
intent(in) :: ia, is, ie
234 integer,
intent(in) :: ja, js, je
235 real(
rp),
intent(in) :: cz ( ka)
236 real(
rp),
intent(in) :: fz (0:ka)
237 real(
rp),
intent(in) :: real_cz ( ka,ia,ja)
238 real(
rp),
intent(in) :: real_fz (0:ka,ia,ja)
239 real(
rp),
intent(in) :: real_phi( ka,ia,ja)
245 log_info(
"ATMOS_REFSTATE_read",*)
'Input reference state profile '
247 if ( atmos_refstate_in_basename /=
'' )
then
251 if ( atmos_refstate_in_check_coordinates )
then
252 call file_cartesc_check_coordinates( fid, atmos=.true. )
256 call file_cartesc_read( fid,
'PRES_ref',
'Z', atmos_refstate1d_pres(:) )
257 call file_cartesc_read( fid,
'TEMP_ref',
'Z', atmos_refstate1d_temp(:) )
258 call file_cartesc_read( fid,
'DENS_ref',
'Z', atmos_refstate1d_dens(:) )
259 call file_cartesc_read( fid,
'POTT_ref',
'Z', atmos_refstate1d_pott(:) )
260 call file_cartesc_read( fid,
'QV_ref',
'Z', atmos_refstate1d_qv(:) )
270 log_error(
"ATMOS_REFSTATE_read",*)
'refstate file is not specified.'
275 cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
287 logical,
save :: first = .true.
290 if ( .not. first )
return
293 if ( atmos_refstate_out_basename /=
'' )
then
296 log_info(
"ATMOS_REFSTATE_write",*)
'Output reference state profile '
299 call file_cartesc_write( atmos_refstate1d_pres(:), atmos_refstate_out_basename, atmos_refstate_out_title, &
300 'PRES_ref',
'Reference profile of pres.',
'Pa',
'Z', atmos_refstate_out_dtype )
301 call file_cartesc_write( atmos_refstate1d_temp(:), atmos_refstate_out_basename, atmos_refstate_out_title, &
302 'TEMP_ref',
'Reference profile of temp.',
'K',
'Z', atmos_refstate_out_dtype )
303 call file_cartesc_write( atmos_refstate1d_dens(:), atmos_refstate_out_basename, atmos_refstate_out_title, &
304 'DENS_ref',
'Reference profile of rho',
'kg/m3',
'Z', atmos_refstate_out_dtype )
305 call file_cartesc_write( atmos_refstate1d_pott(:), atmos_refstate_out_basename, atmos_refstate_out_title, &
306 'POTT_ref',
'Reference profile of theta',
'K',
'Z', atmos_refstate_out_dtype )
307 call file_cartesc_write( atmos_refstate1d_qv(:), atmos_refstate_out_basename, atmos_refstate_out_title, &
308 'QV_ref',
'Reference profile of qv',
'kg/kg',
'Z', atmos_refstate_out_dtype )
311 call file_cartesc_write(
atmos_refstate_pres(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, &
312 'PRES_ref3D',
'Reference profile of pres.',
'Pa',
'ZXY', atmos_refstate_out_dtype )
313 call file_cartesc_write(
atmos_refstate_temp(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, &
314 'TEMP_ref3D',
'Reference profile of temp.',
'K',
'ZXY', atmos_refstate_out_dtype )
315 call file_cartesc_write(
atmos_refstate_dens(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, &
316 'DENS_ref3D',
'Reference profile of rho',
'kg/m3',
'ZXY', atmos_refstate_out_dtype )
317 call file_cartesc_write(
atmos_refstate_pott(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, &
318 'POTT_ref3D',
'Reference profile of theta',
'K',
'ZXY', atmos_refstate_out_dtype )
319 call file_cartesc_write(
atmos_refstate_qv(:,:,:), atmos_refstate_out_basename, atmos_refstate_out_title, &
320 'QV_ref3D',
'Reference profile of qv',
'kg/kg',
'ZXY', atmos_refstate_out_dtype )
329 subroutine atmos_refstate_generate_isa( &
330 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
331 CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
336 profile_isa => atmos_profile_isa
338 hydrostatic_buildrho => atmos_hydrostatic_buildrho
340 saturation_psat_all => atmos_saturation_psat_all, &
341 saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
343 integer,
intent(in) :: ka, ks, ke
344 integer,
intent(in) :: ia, is, ie
345 integer,
intent(in) :: ja, js, je
347 real(
rp),
intent(in) :: cz ( ka)
348 real(
rp),
intent(in) :: fz (0:ka)
349 real(
rp),
intent(in) :: real_cz ( ka,ia,ja)
350 real(
rp),
intent(in) :: real_fz (0:ka,ia,ja)
351 real(
rp),
intent(in) :: real_phi( ka,ia,ja)
372 pott_sfc = atmos_refstate_temp_sfc
375 call profile_isa( ka, ks, ke, &
376 pott_sfc, pres_sfc, &
386 call hydrostatic_buildrho( ka, ks, ke, &
387 pott(:), qv(:), qc(:), &
388 pres_sfc, pott_sfc, qv_sfc, qc_sfc, &
390 dens(:), temp(:), pres(:), &
394 call saturation_psat_all( temp_sfc, psat_sfc )
395 call saturation_dens2qsat_all( ka, ks, ke, &
399 psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc
400 qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp-epsvap) * psat_sfc )
402 qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
406 call hydrostatic_buildrho( ka, ks, ke, &
407 pott(:), qv(:), qc(:), &
408 pres_sfc, pott_sfc, qv_sfc, qc_sfc, &
410 dens(:), temp(:), pres(:), &
413 atmos_refstate1d_pres(:) = pres(:)
414 atmos_refstate1d_temp(:) = temp(:)
415 atmos_refstate1d_dens(:) = dens(:)
416 atmos_refstate1d_pott(:) = pott(:)
417 atmos_refstate1d_qv(:) = qv(:)
420 cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
423 end subroutine atmos_refstate_generate_isa
427 subroutine atmos_refstate_generate_uniform( &
428 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
429 CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
434 hydrostatic_buildrho => atmos_hydrostatic_buildrho
436 saturation_psat_all => atmos_saturation_psat_all, &
437 saturation_dens2qsat_all => atmos_saturation_dens2qsat_all
439 integer,
intent(in) :: ka, ks, ke
440 integer,
intent(in) :: ia, is, ie
441 integer,
intent(in) :: ja, js, je
443 real(
rp),
intent(in) :: cz ( ka)
444 real(
rp),
intent(in) :: fz (0:ka)
445 real(
rp),
intent(in) :: real_cz ( ka,ia,ja)
446 real(
rp),
intent(in) :: real_fz (0:ka,ia,ja)
447 real(
rp),
intent(in) :: real_phi( ka,ia,ja)
469 pott_sfc = atmos_refstate_temp_sfc
474 pott(k) = atmos_refstate_pott_uniform
480 call hydrostatic_buildrho( ka, ks, ke, &
481 pott(:), qv(:), qc(:), &
482 pres_sfc, pott_sfc, qv_sfc, qc_sfc, &
484 dens(:), temp(:), pres(:), &
488 call saturation_psat_all( temp_sfc, psat_sfc )
489 call saturation_dens2qsat_all( ka, ks, ke, &
493 psat_sfc = atmos_refstate_rh * 1.e-2_rp * psat_sfc
494 qv_sfc = epsvap * psat_sfc / ( pres_sfc - (1.0_rp - epsvap) * psat_sfc )
496 qv(k) = atmos_refstate_rh * 1.e-2_rp * qsat(k)
500 call hydrostatic_buildrho( ka, ks, ke, &
501 pott(:), qv(:), qc(:), &
502 pres_sfc, pott_sfc, qv_sfc, qc_sfc, &
504 dens(:), temp(:), pres(:), &
507 atmos_refstate1d_pres(:) = pres(:)
508 atmos_refstate1d_temp(:) = temp(:)
509 atmos_refstate1d_dens(:) = dens(:)
510 atmos_refstate1d_pott(:) = pott(:)
511 atmos_refstate1d_qv(:) = qv(:)
514 cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
517 end subroutine atmos_refstate_generate_uniform
521 subroutine atmos_refstate_generate_zero( KA, IA, JA )
523 integer,
intent(in) :: ka, ia, ja
529 atmos_refstate1d_pres(k) = 0.0_rp
530 atmos_refstate1d_temp(k) = 0.0_rp
531 atmos_refstate1d_dens(k) = 0.0_rp
532 atmos_refstate1d_pott(k) = 0.0_rp
533 atmos_refstate1d_qv(k) = 0.0_rp
549 end subroutine atmos_refstate_generate_zero
554 KA, KS, KE, IA, IS, IE, ISB, IEB, JA, JS, JE, JSB, JEB, &
555 DENS, POTT, TEMP, PRES, QV, &
556 CZ, FZ, FDZ, RCDZ, REAL_CZ, REAL_FZ, REAL_PHI, AREA, &
560 statistics_horizontal_mean
565 integer,
intent(in) :: ka, ks, ke
566 integer,
intent(in) :: ia, is, ie, isb, ieb
567 integer,
intent(in) :: ja, js, je, jsb, jeb
568 real(
rp),
intent(in) :: dens (ka,ia,ja)
569 real(
rp),
intent(in) :: pott (ka,ia,ja)
570 real(
rp),
intent(in) :: temp (ka,ia,ja)
571 real(
rp),
intent(in) :: pres (ka,ia,ja)
572 real(
rp),
intent(in) :: qv (ka,ia,ja)
573 real(
rp),
intent(in) :: cz ( ka)
574 real(
rp),
intent(in) :: fz (0:ka)
575 real(
rp),
intent(in) :: fdz ( ka-1)
576 real(
rp),
intent(in) :: rcdz ( ka)
577 real(
rp),
intent(in) :: real_cz ( ka,ia,ja)
578 real(
rp),
intent(in) :: real_fz (0:ka,ia,ja)
579 real(
rp),
intent(in) :: real_phi( ka,ia,ja)
580 real(
rp),
intent(in) :: area ( ia,ja)
581 real(
dp),
intent(in) :: nowsec
582 logical,
intent(in),
optional :: force
584 real(
rp) :: work(ka,ia,ja)
590 if (
present(force) )
then
596 if ( force_ .or. ( nowsec - last_updated >= atmos_refstate_update_dt ) )
then
598 log_info(
"ATMOS_REFSTATE_update",*)
'update reference state'
601 cz(:), real_cz(:,:,:), temp(:,:,:), &
603 call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
604 work(:,:,:), area(:,:), &
605 atmos_refstate1d_temp(:) )
608 cz(:), real_cz(:,:,:), pres(:,:,:), &
610 call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
611 work(:,:,:), area(:,:), &
612 atmos_refstate1d_pres(:) )
615 cz(:), real_cz(:,:,:), dens(:,:,:), &
617 call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
618 work(:,:,:), area(:,:), &
619 atmos_refstate1d_dens(:) )
622 cz(:), real_cz(:,:,:), pott(:,:,:), &
624 call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
625 work(:,:,:), area(:,:), &
626 atmos_refstate1d_pott(:) )
629 cz(:), real_cz(:,:,:), qv(:,:,:), &
631 call statistics_horizontal_mean( ka, ks, ke, ia, is, ie, ja, js, je, &
632 work(:,:,:), area(:,:), &
633 atmos_refstate1d_qv(:) )
636 if( atmos_refstate1d_dens(k) <= 0.0_rp ) atmos_refstate1d_dens(k) = atmos_refstate1d_dens(k+1)
637 if( atmos_refstate1d_temp(k) <= 0.0_rp ) atmos_refstate1d_temp(k) = atmos_refstate1d_temp(k+1)
638 if( atmos_refstate1d_pres(k) <= 0.0_rp ) atmos_refstate1d_pres(k) = atmos_refstate1d_pres(k+1)
639 if( atmos_refstate1d_pott(k) <= 0.0_rp ) atmos_refstate1d_pott(k) = atmos_refstate1d_pott(k+1)
640 if( atmos_refstate1d_qv(k) <= 0.0_rp ) atmos_refstate1d_qv(k) = atmos_refstate1d_qv(k+1)
647 cz(:), fz(:), real_cz(:,:,:), real_fz(:,:,:), real_phi(:,:,:) )
649 last_updated = nowsec
653 log_info(
"ATMOS_REFSTATE_update",*)
'Generated reference state of atmosphere:'
654 log_info_cont(*)
'============================================================================='
655 log_info_cont(*)
' z*-coord.: pressure: temperature: density: pot.temp.: water vapor'
657 log_info_cont(
'(6F13.5)') cz(k), &
658 atmos_refstate1d_pres(k), &
659 atmos_refstate1d_temp(k), &
660 atmos_refstate1d_dens(k), &
661 atmos_refstate1d_pott(k), &
662 atmos_refstate1d_qv(k)
664 log_info_cont(*)
'============================================================================='
677 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
678 CZ, FZ, REAL_CZ, REAL_FZ, REAL_PHI )
690 hydrostatic_buildrho_atmos => atmos_hydrostatic_buildrho_atmos, &
691 hydrostatic_buildrho_atmos_rev => atmos_hydrostatic_buildrho_atmos_rev
693 integer,
intent(in) :: KA, KS, KE
694 integer,
intent(in) :: IA, IS, IE
695 integer,
intent(in) :: JA, JS, JE
697 real(RP),
intent(in) :: CZ ( KA)
698 real(RP),
intent(in) :: FZ (0:KA)
699 real(RP),
intent(in) :: REAL_CZ ( KA,IA,JA)
700 real(RP),
intent(in) :: REAL_FZ (0:KA,IA,JA)
701 real(RP),
intent(in) :: REAL_PHI( KA,IA,JA)
703 real(RP) :: dens(KA,IA,JA)
704 real(RP) :: temp(KA,IA,JA)
705 real(RP) :: pres(KA,IA,JA)
706 real(RP) :: pott(KA,IA,JA)
707 real(RP) :: qv (KA,IA,JA)
708 real(RP) :: qc (KA,IA,JA)
709 real(RP) :: dz (KA,IA,JA)
711 real(RP) :: dens_toa_1D
712 real(RP) :: temp_toa_1D
713 real(RP) :: pres_toa_1D
717 real(RP) :: work(KA,IA,JA)
727 work(:,i,j) = atmos_refstate1d_pott(:)
732 real_cz(:,:,:), cz(:), &
739 work(:,i,j) = atmos_refstate1d_qv(:)
744 real_cz(:,:,:), cz(:), &
751 dz_1d = fz(ke) - cz(ke)
753 call hydrostatic_buildrho_atmos( atmos_refstate1d_pott(ke), &
754 atmos_refstate1d_qv(ke), &
756 atmos_refstate1d_dens(ke), &
757 atmos_refstate1d_pott(ke), &
758 atmos_refstate1d_qv(ke), &
769 dz(ks-1,i,j) = real_cz(ks,i,j) - real_fz(ks-1,i,j)
771 dz(k,i,j) = real_cz(k+1,i,j) - real_cz(k,i,j)
773 dz(ke,i,j) = real_fz(ke,i,j) - real_cz(ke,i,j)
779 dens(ke+1,i,j) = dens_toa_1d
780 temp(ke+1,i,j) = temp_toa_1d
781 pres(ke+1,i,j) = pres_toa_1d
782 pott(ke+1,i,j) = pott(ke,i,j)
783 qv(ke+1,i,j) = qv(ke,i,j)
789 pott(ks-1,i,j) = pott(ks,i,j)
790 qv(ks-1,i,j) = qv(ks,i,j)
796 call hydrostatic_buildrho_atmos_rev( ia, is, ie, ja, js, je, &
797 pott(ke,:,:), qv(ke,:,:), qc(ke,:,:), &
798 dens(ke+1,:,:), pott(ke+1,:,:), qv(ke+1,:,:), qc(ke+1,:,:), &
800 dens(ke ,:,:), temp(ke ,:,:), pres(ke ,:,:) )
802 call hydrostatic_buildrho_atmos_rev( ka, ks, ke, ia, is, ie, ja, js, je, &
803 pott(:,:,:), qv(:,:,:), qc(:,:,:), &
806 temp(:,:,:), pres(:,:,:) )
872 integer,
intent(in) :: KA, KS, KE
874 real(RP),
intent(in) :: FDZ (KA-1)
875 real(RP),
intent(in) :: RCDZ(KA)
877 real(RP),
intent(inout) :: phi(KA)
880 real(RP) :: correction(KA)
883 integer,
parameter :: iter_max = 100
884 real(RP) :: sig0, sig1, zerosw
892 correction(ks-1:ks+1) = 0.0_rp
893 correction(ke-1:ke+1) = 0.0_rp
895 fact(ks-1:ks+1) = 0.0_rp
896 fact(ke-1:ke+1) = 0.0_rp
898 do iter = 1, iter_max
902 dev(k) = phi(k) - ( fdz(k-1)*phi(k+1) + fdz(k)*phi(k-1) ) / ( fdz(k) + fdz(k-1) )
906 sig0 = dev(k) * dev(k-1)
907 sig1 = dev(k) * dev(k+1)
908 if ( sig0 < -eps .and. sig1 < -eps )
then
909 correction(k) = dev(k) &
910 / ( 2.0_rp*rcdz(k) + ( fdz(k-1)*rcdz(k+1) + fdz(k)*rcdz(k-1) ) / ( fdz(k) + fdz(k-1) ) )
913 correction(k) = 0.0_rp
917 sig1 = dev(ks+1) * dev(ks+2)
918 if ( sig1 < -eps )
then
919 correction(ks+1) = dev(ks+1) &
920 / ( 2.0_rp*rcdz(ks+1) + (fdz(ks)*rcdz(ks+2)+fdz(ks+1)*rcdz(ks))/(fdz(ks+1)+fdz(ks)) )
923 correction(ks+1) = 0.0_rp
926 sig0 = dev(ke-1) * dev(ke-2)
927 if ( sig0 < -eps )
then
928 correction(ke-1) = dev(ke-1) &
929 / ( 2.0_rp*rcdz(ke-1) + (fdz(ke-2)*rcdz(ke)+fdz(ke-1)*rcdz(ke-2))/(fdz(ke-1)+fdz(ke-2)) )
932 correction(ke-1) = 0.0_rp
935 if ( .NOT. updated )
exit
938 zerosw = 0.5_rp - sign( 0.5_rp, abs(correction(k))-eps )
939 fact(k) = correction(k) / ( correction(k) - correction(k+1) - correction(k-1) + zerosw )
943 phi(k) = phi(k) + ( correction(k+1) * fact(k+1) &
944 - correction(k ) * fact(k ) * 2.0_rp &
945 + correction(k-1) * fact(k-1) ) * rcdz(k)
948 if ( iter == iter_max )
then
949 log_info(
"ATMOS_REFSTATE_smoothing",*)
"iteration not converged!", phi