Changeset 820 for codes/icosagcm/devel/src/dissip/nudging_mod.f90
- Timestamp:
- 04/02/19 18:10:35 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dissip/nudging_mod.f90
r818 r820 1 1 MODULE nudging_mod 2 2 USE icosa, ONLY : rstd 3 USE grid_param, ONLY : llm 3 USE grid_param, ONLY : llm, nqdyn 4 USE omp_para, ONLY : ll_begin, ll_end 4 5 USE domain_mod, ONLY : ndomain, assigned_domain 5 6 USE dimensions, ONLY : swap_dimensions, u_right, u_lup, u_ldown … … 12 13 ! nudging will be active outside a disc of radius 'radius' centered at 'center_lon', 'center lat'. 13 14 REAL(rstd) :: center_lon, center_lat, nudging_radius, time 15 !$OMP THREADPRIVATE(center_lon, center_lat, nudging_radius, time) 14 16 15 TYPE(t_field),POINTER :: f_relax_coef_e(:), f_target_ue(:) 17 TYPE(t_field),POINTER :: f_relax_coef_e(:), f_target_ue(:), & 18 f_relax_coef_i(:), f_target_theta_rhodz(:) 16 19 17 20 PUBLIC :: init_guided, guided … … 19 22 CONTAINS 20 23 21 SUBROUTINE init_guided(f_u e)24 SUBROUTINE init_guided(f_u,f_theta_rhodz) 22 25 USE getin_mod, ONLY : getin 23 26 USE math_const, ONLY : pi 24 27 USE earth_const, ONLY : scale_factor 25 TYPE(t_field),POINTER :: f_ue(:)! initial condition 28 TYPE(t_field),POINTER :: f_u(:)! initial condition 29 TYPE(t_field),POINTER :: f_theta_rhodz(:)! initial condition 26 30 REAL(rstd), POINTER :: ue(:,:), target_ue(:,:), coef_e(:) 31 REAL(rstd), POINTER :: theta_rhodz(:,:,:), target_theta_rhodz(:,:,:), coef_i(:) 27 32 INTEGER :: ind 28 33 ! read DEF keys describing how to relax … … 39 44 ! we should check that radius>0 40 45 41 CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e') 42 CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e') 46 !SELECT CASE(TRIM(nudg_name)) 47 ! CASE ('f_u') 48 CALL allocate_field(f_relax_coef_e, field_u, type_real, name='nudging_coef_e') 49 CALL allocate_field(f_target_ue, field_u, type_real, llm, name='nudging_target_e') 50 ! CASE ('f_theta_rhodz') 51 CALL allocate_field(f_relax_coef_i, field_t, type_real, name='nudging_coef_i') 52 CALL allocate_field(f_target_theta_rhodz, field_t, type_real, llm,nqdyn, name='nudging_target_theta') 53 ! CASE DEFAULT 54 !END SELECT 43 55 44 56 ! compute relax_coef and target_ue … … 46 58 center_lat = center_lat * pi/180. 47 59 DO ind = 1 , ndomain 48 IF (.NOT. assigned_domain(ind)) CYCLE 49 CALL swap_dimensions(ind) 50 CALL swap_geometry(ind) 51 coef_e = f_relax_coef_e(ind) 52 target_ue = f_target_ue(ind) 53 ue = f_ue(ind) 54 CALL compute_relax_coef(coef_e) 55 CALL compute_target(ue, target_ue) 60 IF (.NOT. assigned_domain(ind)) CYCLE 61 CALL swap_dimensions(ind) 62 CALL swap_geometry(ind) 63 coef_e = f_relax_coef_e(ind) 64 coef_i = f_relax_coef_i(ind) 65 target_ue = f_target_ue(ind) 66 ue = f_u(ind) 67 target_theta_rhodz = f_target_theta_rhodz(ind) 68 theta_rhodz = f_theta_rhodz(ind) 69 CALL compute_relax_coef(coef_e, coef_i) 70 CALL compute_target_u(ue, target_ue) 71 CALL compute_target_center(theta_rhodz, target_theta_rhodz) 56 72 END DO 57 73 END SUBROUTINE init_guided 74 75 !----------------------------- Compute relaxation coefficients ------------------------------ 58 76 59 SUBROUTINE compute_relax_coef(coef_e) 60 USE omp_para, ONLY : ll_begin, ll_end 77 SUBROUTINE compute_relax_coef(coef_e, coef_i) 61 78 USE geometry, ONLY : lon_e, lat_e, lon_i, lat_i 62 REAL(rstd), INTENT(OUT) :: coef_e(iim*3*jjm) 79 REAL(rstd), INTENT(OUT) :: coef_e(iim*3*jjm), coef_i(iim*jjm) 63 80 INTEGER :: l, ij 64 coef_e(:)=1.65 81 DO ij=ij_begin_ext, ij_end_ext 66 82 coef_e(ij+u_right) = relax_coef(lon_e(ij+u_right), lat_e(ij+u_right) ) 67 83 coef_e(ij+u_lup) = relax_coef(lon_e(ij+u_lup), lat_e(ij+u_lup) ) 68 84 coef_e(ij+u_ldown) = relax_coef(lon_e(ij+u_ldown), lat_e(ij+u_ldown) ) 85 coef_i(ij) = relax_coef(lon_i(ij), lat_i(ij) ) 69 86 END DO 70 87 END SUBROUTINE compute_relax_coef 71 88 72 89 FUNCTION relax_coef(lon,lat) 73 90 USE spherical_geom_mod, ONLY : dist_lonlat … … 77 94 CALL dist_lonlat(lon, lat, center_lon, center_lat, dist) 78 95 c = tanh((1.-radius*dist/nudging_radius)*20.) ! 1 inside circle, -1 outside 79 !c = tanh((1.-radius*dist/nudging_radius)*10.) ! 1 inside circle, -1 outside80 !c = tanh((1.-radius*dist/nudging_radius)*5.) ! 1 inside circle, -1 outside81 96 relax_coef = .5*(1.+c) ! rescale to [0,1] range ; ! 1 inside circle, 0 outside 82 !relax_coef = 1. ! rescale to [0,1] range ; ! 1 inside circle, 0 outside83 PRINT *, 'c,lat,dist,radius,nudging_radius,relax_coef',c,lat,dist,radius,nudging_radius,relax_coef84 97 END FUNCTION relax_coef 85 98 86 SUBROUTINE compute_target(ue, target_ue) 87 USE omp_para, ONLY : ll_begin, ll_end 99 !----------------------------- Copy initial condition as target ------------------------------ 100 101 SUBROUTINE compute_target_u(ue, target_ue) 88 102 REAL(rstd), INTENT(OUT) :: target_ue(iim*3*jjm,llm) 89 103 REAL(rstd), INTENT(IN) :: ue(iim*3*jjm,llm) … … 96 110 END DO 97 111 END DO 98 END SUBROUTINE compute_target 112 END SUBROUTINE compute_target_u 113 114 SUBROUTINE compute_target_center(theta_rhodz, target_theta_rhodz) 115 REAL(rstd), INTENT(OUT) :: target_theta_rhodz(iim*jjm,llm,nqdyn) 116 REAL(rstd), INTENT(IN) :: theta_rhodz(iim*jjm,llm,nqdyn) 117 INTEGER :: l, ij, iq 118 DO iq=1, nqdyn 119 DO l = ll_begin, ll_end 120 DO ij=ij_begin_ext, ij_end_ext 121 target_theta_rhodz(ij,l,iq)=theta_rhodz(ij,l,iq) 122 END DO 123 END DO 124 END DO 125 END SUBROUTINE compute_target_center 126 127 !----------------------------- Relax towards target ------------------------------ 99 128 100 129 SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q) … … 105 134 TYPE(t_field),POINTER :: f_u(:) 106 135 TYPE(t_field),POINTER :: f_q(:) 107 REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:) 136 REAL(rstd), POINTER :: target_ue(:,:), ue(:,:), coef_e(:) 137 REAL(rstd), POINTER :: target_theta_rhodz(:,:,:), theta_rhodz(:,:,:), coef_i(:) 108 138 INTEGER :: ind 109 139 … … 113 143 CALL swap_geometry(ind) 114 144 coef_e = f_relax_coef_e(ind) 115 PRINT *, 'max coeff', maxval(coef_e)116 ! PRINT *, 'coeff', coef_e117 145 target_ue = f_target_ue(ind) 118 146 ue = f_u(ind) 119 CALL compute_guided(coef_e, target_ue, ue) 147 CALL compute_guided_u(coef_e, target_ue, ue) 148 coef_i = f_relax_coef_i(ind) 149 target_theta_rhodz = f_target_theta_rhodz(ind) 150 theta_rhodz = f_theta_rhodz(ind) 151 CALL compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz) 120 152 END DO 121 153 122 154 END SUBROUTINE guided 123 155 124 SUBROUTINE compute_guided(coef_e, target_ue, ue) 125 USE omp_para, ONLY : ll_begin, ll_end 156 SUBROUTINE compute_guided_u(coef_e, target_ue, ue) 126 157 REAL(rstd), INTENT(IN) :: coef_e(iim*3*jjm) 127 158 REAL(rstd), INTENT(IN) :: target_ue(iim*3*jjm,llm) … … 138 169 END DO 139 170 END DO 140 END SUBROUTINE compute_guided 171 END SUBROUTINE compute_guided_u 172 173 SUBROUTINE compute_guided_center(coef_i, target_theta_rhodz, theta_rhodz) 174 REAL(rstd), INTENT(IN) :: coef_i(iim*jjm) 175 REAL(rstd), INTENT(IN) :: target_theta_rhodz(iim*jjm,llm,nqdyn) 176 REAL(rstd), INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn) 177 INTEGER :: l, ij, iq 178 DO iq=1, nqdyn 179 DO l = ll_begin, ll_end 180 DO ij=ij_begin_ext, ij_end_ext 181 theta_rhodz(ij,l,iq) = theta_rhodz(ij,l,iq)*coef_i(ij) + & 182 target_theta_rhodz(ij,l,iq)*(1.-coef_i(ij)) 183 END DO 184 END DO 185 END DO 186 END SUBROUTINE compute_guided_center 141 187 142 188 END MODULE nudging_mod
Note: See TracChangeset
for help on using the changeset viewer.