Changeset 856 for codes/icosagcm/devel/src
- Timestamp:
- 05/06/19 01:49:12 (5 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/base/dimensions.f90
r533 r856 80 80 81 81 TYPE(t_domain),POINTER :: d 82 82 83 IF(.NOT. swap_needed) RETURN 84 83 85 d=>domain(ind) 84 86 -
codes/icosagcm/devel/src/dissip/dissip_gcm.f90
r845 r856 123 123 CALL getin("niterdivgrad",niterdivgrad) 124 124 125 CALL dissip_constants 126 CALL dissip_profile 125 IF(grid_type == grid_ico) THEN 126 CALL dissip_constants 127 CALL dissip_profile 128 END IF 129 127 130 END SUBROUTINE init_dissip 128 131 -
codes/icosagcm/devel/src/icosa_init.f90
r846 r856 42 42 !$OMP PARALLEL 43 43 CALL switch_omp_no_distrib_level 44 IF(grid_type == grid_ico) CALL compute_geometry 44 IF(grid_type == grid_ico) THEN 45 CALL compute_geometry 46 ELSE 47 swap_needed=.FALSE. ! do this inside an OMP PARALLEL loop to affect all threads 48 END IF 45 49 CALL check_total_area 46 50 -
codes/icosagcm/devel/src/initial/etat0_collocated.f90
r848 r856 1 1 MODULE etat0_collocated_mod 2 2 USE icosa 3 USE grid_param 3 4 USE omp_para, ONLY : is_omp_level_master 4 5 USE caldyn_vars_mod, ONLY : hydrostatic … … 6 7 PRIVATE 7 8 8 LOGICAL :: autoinit_mass, autoinit_NH9 9 CHARACTER(len=255),SAVE :: etat0_type 10 !$OMP THREADPRIVATE( autoinit_mass, autoinit_NH,etat0_type)10 !$OMP THREADPRIVATE(etat0_type) 11 11 12 12 PUBLIC :: etat0_type, etat0_collocated … … 36 36 TYPE(t_field),POINTER :: f_q(:) 37 37 38 TYPE(t_field),POINTER,SAVE :: f_temp(:) 38 TYPE(t_field),POINTER,SAVE :: f_temp(:) ! SAVE => all threads see the same f_temp 39 39 REAL(rstd),POINTER :: ps(:) 40 40 REAL(rstd),POINTER :: mass(:,:) … … 64 64 q=f_q(ind) 65 65 66 IF( TRIM(etat0_type)=='williamson91.6' ) THEN 67 CALL compute_etat0_collocated_hex(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q) 68 ELSE 69 CALL compute_etat0_collocated_hex(ps,mass, phis, temp, u, geopot, W, q) 70 ENDIF 71 72 IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) 73 66 SELECT CASE(grid_type) 67 CASE(grid_ico) 68 CALL compute_etat0_collocated_hex(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 69 CASE(grid_unst) 70 CALL compute_etat0_collocated_unst(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q) 71 CASE DEFAULT 72 STOP 'Unexpected value of grid_type encountered in etat0_collocated.' 73 END SELECT 74 74 75 ENDDO 75 76 … … 78 79 END SUBROUTINE etat0_collocated 79 80 80 SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) 81 USE icosa 82 USE pression_mod 83 USE exner_mod 84 USE omp_para 85 REAL(rstd),INTENT(IN) :: ps(iim*jjm) 86 REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm) 87 REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) 88 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 89 INTEGER,INTENT(IN) :: offset 90 91 REAL(rstd) :: p(iim*jjm,llm+1) 92 REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta 93 INTEGER :: i,j,ij,l 94 95 cppd=cpp 96 Rd=kappa*cppd 97 98 CALL compute_pression(ps,p,offset) 99 ! flush p 100 DO l = ll_begin, ll_end 101 DO j=jj_begin-offset,jj_end+offset 102 DO i=ii_begin-offset,ii_end+offset 103 ij=(j-1)*iim+i 104 mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass 105 p_ij = .5*(p(ij,l)+p(ij,l+1)) ! pressure at full level 106 SELECT CASE(caldyn_thermo) 107 CASE(thermo_theta) 108 theta = temp(ij,l)*(p_ij/preff)**(-kappa) 109 theta_rhodz(ij,l) = mass * theta 110 CASE(thermo_entropy) 111 nu = log(p_ij/preff) 112 chi = log(temp(ij,l)/Treff) 113 entropy = cppd*chi-Rd*nu 114 theta_rhodz(ij,l) = mass * entropy 115 ! CASE(thermo_moist) 116 ! q_ij=q(ij,l,1) 117 ! r_ij=1.-q_ij 118 ! mass=mass*(1-q_ij) ! dry mass 119 ! nu = log(p_ij/preff) 120 ! chi = log(temp(ij,l)/Treff) 121 ! entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) 122 ! theta_rhodz(ij,l) = mass * entropy 123 CASE DEFAULT 124 STOP 125 END SELECT 126 ENDDO 127 ENDDO 128 ENDDO 129 END SUBROUTINE compute_temperature2entropy 130 131 SUBROUTINE compute_etat0_collocated_hex(ps,mass,phis,temp_i,u, geopot,W, q) 132 USE wind_mod 133 USE disvert_mod 81 SUBROUTINE compute_etat0_collocated_hex(phis,ps,mass,theta_rhodz,u, geopot,W, q) 134 82 REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 135 83 REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm) 84 REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 136 85 REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 137 REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)138 86 REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 139 87 REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) … … 141 89 REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 142 90 143 REAL(rstd) :: ulon_i(iim*jjm,llm)144 REAL(rstd) :: ulat_i(iim*jjm,llm)145 146 91 REAL(rstd) :: ps_e(3*iim*jjm) 92 REAL(rstd) :: phis_e(3*iim*jjm) 147 93 REAL(rstd) :: mass_e(3*iim*jjm,llm) 148 REAL(rstd) :: phis_e(3*iim*jjm)149 REAL(rstd) :: temp_e(3*iim*jjm,llm)150 94 REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) 151 REAL(rstd) :: ulon_e(3*iim*jjm,llm)152 REAL(rstd) :: ulat_e(3*iim*jjm,llm)153 95 REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot) 154 96 155 INTEGER :: l,i,j,ij156 REAL :: p_ik, v_ik, mass_ik157 158 ! For NH geopotential and vertical momentum must be initialized.159 ! Unless autoinit_NH is set to .FALSE. , they will be initialized160 ! to hydrostatic geopotential and zero161 autoinit_mass = .TRUE.162 autoinit_NH = .NOT. hydrostatic163 97 w(:,:) = 0 164 165 CALL compute_etat0_collocated(iim*jjm , lon_i, lat_i, phis, ps, mass, temp_i, ulon_i, ulat_i, geopot, q) 166 CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 167 168 IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 169 IF(autoinit_NH) THEN 170 geopot(:,1) = phis(:) ! surface geopotential 171 DO l = 1, llm 172 DO ij=1,iim*jjm 173 ! hybrid pressure coordinate 174 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 175 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 176 ! v=R.T/p, R=kappa*cpp 177 v_ik = kappa*cpp*temp_i(ij,l)/p_ik 178 geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 179 END DO 180 END DO 181 END IF 182 183 CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u) 98 CALL compute_etat0_collocated(iim*jjm , lon_i, lat_i, phis, ps, mass, theta_rhodz, geopot, q) 99 CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u) 184 100 185 101 END SUBROUTINE compute_etat0_collocated_hex 186 102 187 SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, temp, ulon, ulat, geopot, q) 103 SUBROUTINE compute_etat0_collocated_unst(phis,ps,mass,theta_rhodz,u, geopot,W, q) 104 REAL(rstd),INTENT(INOUT) :: ps(primal_num) 105 REAL(rstd),INTENT(INOUT) :: mass(llm,primal_num) 106 REAL(rstd),INTENT(OUT) :: theta_rhodz(llm,primal_num) 107 REAL(rstd),INTENT(OUT) :: phis(primal_num) 108 REAL(rstd),INTENT(OUT) :: u(llm, edge_num) 109 REAL(rstd),INTENT(OUT) :: W(llm+1, primal_num) 110 REAL(rstd),INTENT(OUT) :: geopot(llm+1, primal_num) 111 REAL(rstd),INTENT(OUT) :: q(llm, primal_num, nqtot) 112 ! The 2D/3D arrays below have the shape expected by compute_etat0_collocated 113 REAL(rstd) :: ps_e(edge_num) 114 REAL(rstd) :: phis_e(edge_num) 115 REAL(rstd) :: u_e(edge_num, llm) 116 REAL(rstd) :: mass_i(primal_num, llm), theta_rhodz_i(primal_num, llm), mass_e(edge_num, llm) 117 REAL(rstd) :: geopot_i(edge_num, llm+1), geopot_e(edge_num, llm+1) 118 REAL(rstd) :: q_i(primal_num, llm, nqtot), q_e(edge_num, llm, nqtot) 119 INTEGER :: iq 120 121 w(:,:) = 0 122 CALL compute_etat0_collocated(primal_num , lon_i, lat_i, phis, ps, mass_i, theta_rhodz_i, geopot_i, q_i) 123 CALL compute_etat0_collocated(edge_num, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u_e) 124 mass = TRANSPOSE(mass_i) 125 theta_rhodz = TRANSPOSE(theta_rhodz_i) 126 geopot = TRANSPOSE(geopot_i) 127 u = TRANSPOSE(u_e) 128 DO iq=1,nqtot 129 q(:,:,iq) = TRANSPOSE(q_i(:,:,iq)) 130 END DO 131 END SUBROUTINE compute_etat0_collocated_unst 132 133 SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, theta_rhodz, geopot, q, ep, uperp) 188 134 USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0 189 135 USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 … … 199 145 USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0 200 146 USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0 147 USE disvert_mod, ONLY : ptop, ap, bp, mass_ak, mass_bk, mass_dak, mass_dbk 201 148 INTEGER :: ngrid 202 149 REAL(rstd),INTENT(IN) :: lon(ngrid), lat(ngrid) 203 150 REAL(rstd),INTENT(INOUT) :: ps(ngrid) 204 151 REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm) 152 REAL(rstd),INTENT(OUT) :: theta_rhodz(ngrid,llm) 205 153 REAL(rstd),INTENT(OUT) :: phis(ngrid) 206 REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)207 REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)208 REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)209 154 REAL(rstd),INTENT(OUT) :: geopot(ngrid,llm+1) 210 155 REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 156 REAL(rstd),INTENT(OUT), OPTIONAL :: ep(ngrid,3), uperp(ngrid,llm) 157 158 REAL(rstd) :: ulon(ngrid,llm) 159 REAL(rstd) :: ulat(ngrid,llm) 160 REAL(rstd) :: temp(ngrid,llm) 161 162 LOGICAL :: autoinit_mass, autoinit_NH 163 INTEGER :: ij,l 164 REAL(rstd) :: clat, slat, clon, slon, u3d(3), p_ik, mass_ik, v_ik 165 REAL(rstd) :: q_ik, r_ik, chi, entropy, theta, log_p_preff 166 167 ! For NH, geopotential and vertical momentum must be initialized. 168 ! Unless autoinit_NH is set to .FALSE. , they will be initialized 169 ! to hydrostatic geopotential and zero 170 autoinit_mass = .TRUE. 171 autoinit_NH = .NOT. hydrostatic 211 172 212 173 SELECT CASE (TRIM(etat0_type)) … … 232 193 ! autoinit_NH = .FALSE. ! compute_bubble initializes geopot 233 194 CASE('williamson91.6') 234 CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), t emp(:,1), ulon(:,1), ulat(:,1))195 CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), theta_rhodz(:,1), ulon(:,1), ulat(:,1)) 235 196 autoinit_mass = .FALSE. ! do not overwrite mass 236 197 CASE('dcmip2016_baroclinic_wave') … … 242 203 END SELECT 243 204 205 IF(PRESENT(uperp)) THEN 206 DO l = 1, llm 207 DO ij=1,ngrid 208 clat = COS(lat(ij)) 209 slat = SIN(lat(ij)) 210 clon = COS(lon(ij)) 211 slon = SIN(lon(ij)) 212 u3d(1) = -clon*slat*ulat(ij,l) - slon*ulon(ij,l) 213 u3d(2) = -slon*slat*ulat(ij,l) + clon*ulon(ij,l) 214 u3d(3) = clat*ulat(ij,l) 215 uperp(ij,l) = SUM(u3d*ep(ij,:)) 216 END DO 217 END DO 218 END IF 219 220 IF(autoinit_mass) THEN 221 DO l = 1, llm 222 DO ij=1,ngrid 223 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 224 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 225 mass(ij,l) = mass_ik 226 SELECT CASE(caldyn_thermo) 227 CASE(thermo_theta) 228 theta = temp(ij,l)*(p_ik/preff)**(-kappa) 229 theta_rhodz(ij,l) = mass_ik * theta 230 CASE(thermo_entropy) 231 log_p_preff = log(p_ik/preff) 232 chi = log(temp(ij,l)/Treff) 233 entropy = cpp*chi-Rd*log_p_preff 234 theta_rhodz(ij,l) = mass_ik * entropy 235 CASE(thermo_moist) 236 q_ik=q(ij,l,1) 237 r_ik=1.-q_ik 238 mass_ik = mass_ik*(1.-q_ik) ! dry mass 239 log_p_preff = log(p_ik/preff) 240 chi = log(temp(ij,l)/Treff) 241 entropy = r_ik*(cpp*chi-Rd*nu) + q_ik*(cppv*chi-Rv*log_p_preff) 242 theta_rhodz(ij,l) = mass_ik * entropy 243 CASE DEFAULT 244 STOP 245 END SELECT 246 END DO 247 END DO 248 END IF 249 250 IF(autoinit_NH) THEN 251 geopot(:,1) = phis(:) ! surface geopotential 252 DO l = 1, llm 253 DO ij=1,ngrid 254 ! hybrid pressure coordinate 255 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 256 mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 257 ! v=R.T/p, R=kappa*cpp 258 v_ik = kappa*cpp*temp(ij,l)/p_ik 259 geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 260 END DO 261 END DO 262 END IF 263 244 264 END SUBROUTINE compute_etat0_collocated 245 265 -
codes/icosagcm/devel/src/output/output_field.f90
r714 r856 20 20 SUBROUTINE output_field_init 21 21 USE getin_mod 22 USE grid_param 22 23 IMPLICIT NONE 23 24 … … 35 36 36 37 IF (xios_output) THEN 37 CALL xios_init_write_field38 IF(grid_type == grid_ico) CALL xios_init_write_field 38 39 ENDIF 39 40 END SUBROUTINE output_field_init -
codes/icosagcm/devel/src/parallel/domain.f90
r821 r856 65 65 INTEGER,SAVE :: ndomain 66 66 INTEGER,SAVE :: ndomain_glo 67 68 LOGICAL :: swap_needed=.TRUE. ! .FALSE. if a thread always computes on the same domain 69 !$OMP THREADPRIVATE(swap_needed) 67 70 68 71 TYPE(t_domain),SAVE,ALLOCATABLE,TARGET :: domain(:) -
codes/icosagcm/devel/src/sphere/geometry.f90
r602 r856 142 142 SUBROUTINE swap_geometry(ind) 143 143 USE field_mod 144 USE domain_mod, ONLY : swap_needed 144 145 IMPLICIT NONE 145 146 INTEGER,INTENT(IN) :: ind 147 148 IF(.NOT. swap_needed) RETURN 149 146 150 !!$OMP MASTER 147 151 Ai=geom%Ai(ind) -
codes/icosagcm/devel/src/transport/advect_tracer.f90
r592 r856 27 27 SUBROUTINE init_advect_tracer 28 28 USE omp_para 29 USE grid_param 29 30 REAL(rstd),POINTER :: tangent(:,:) 30 31 REAL(rstd),POINTER :: normal(:,:) … … 42 43 CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq') 43 44 44 DO ind=1,ndomain 45 IF (.NOT. assigned_domain(ind)) CYCLE 46 CALL swap_dimensions(ind) 47 CALL swap_geometry(ind) 48 normal=f_normal(ind) 49 tangent=f_tangent(ind) 50 sqrt_leng=f_sqrt_leng(ind) 51 IF (is_omp_level_master) CALL init_advect(normal,tangent,sqrt_leng) 52 END DO 45 IF(grid_type == grid_ico) THEN 46 DO ind=1,ndomain 47 IF (.NOT. assigned_domain(ind)) CYCLE 48 CALL swap_dimensions(ind) 49 CALL swap_geometry(ind) 50 normal=f_normal(ind) 51 tangent=f_tangent(ind) 52 sqrt_leng=f_sqrt_leng(ind) 53 IF (is_omp_level_master) CALL init_advect(normal,tangent,sqrt_leng) 54 END DO 55 END IF 53 56 54 57 END SUBROUTINE init_advect_tracer -
codes/icosagcm/devel/src/unstructured/init_unstructured.f90
r830 r856 96 96 USE ioipsl 97 97 USE field_mod 98 USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e 98 99 IMPLICIT NONE 99 100 100 101 !!-------------Declare local variables------------------- 101 CHARACTER(LEN=*),PARAMETER :: f=" mesh_information.nc"102 INTEGER :: id_nc, ierr, status, descriptionLength 102 CHARACTER(LEN=*),PARAMETER :: f="input/mesh_information.nc" 103 INTEGER :: id_nc, ierr, status, descriptionLength, ij 103 104 CHARACTER(LEN= 80) :: description 104 105 REAL(rstd), ALLOCATABLE :: angle_e(:) 106 REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat 107 elon(3), elat(3) ! lon/lat unit vectors 105 108 print *,"------------------ READING FILE " , f, "----------------------- " 106 109 !open and read the input file … … 167 170 ALLOCATE(dual_vertex, source = Idata_read2) 168 171 172 CALL read_from_gridfile(id_nc, 'float', 'lon_i') 173 ALLOCATE(lon_i, source = Ddata_read1) 174 CALL read_from_gridfile(id_nc, 'float', 'lat_i') 175 ALLOCATE(lat_i, source = Ddata_read1) 176 CALL read_from_gridfile(id_nc, 'float', 'lon_e') 177 ALLOCATE(lon_e, source = Ddata_read1) 178 CALL read_from_gridfile(id_nc, 'float', 'lat_e') 179 ALLOCATE(lat_e, source = Ddata_read1) 180 CALL read_from_gridfile(id_nc, 'float', 'angle_e') 181 ALLOCATE(angle_e, source = Ddata_read1) 182 169 183 CALL free_data_read ! free buffers after reading all data from grid file 170 184 … … 176 190 max_trisk_deg = SIZE(trisk,1) 177 191 192 ! now post-process some of the data we just read 193 ALLOCATE(ep_e(edge_num,3)) 194 DO ij=1,edge_num 195 clon = COS(lon_e(ij)) 196 slon = SIN(lon_e(ij)) 197 clat = COS(lat_e(ij)) 198 slat = SIN(lat_e(ij)) 199 ! x,y,z = clat*clon, clat*slon, slat 200 elon = (/ -slon, clon, 0. /) 201 elat = (/ -slat*clon, -slat*slon, clat /) 202 ep_e(ij,:) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat 203 END DO 204 205 DEALLOCATE(angle_e) 178 206 END SUBROUTINE read_dump_partition 179 207
Note: See TracChangeset
for help on using the changeset viewer.