MODULE etat0_mod USE icosa USE omp_para IMPLICIT NONE PRIVATE CHARACTER(len=255),SAVE :: etat0_type !$OMP THREADPRIVATE(etat0_type) REAL(rstd) :: etat0_temp PUBLIC :: etat0, init_etat0, etat0_type ! Important notes for OpenMP ! When etat0 is called, vertical OpenMP parallelism is deactivated. ! Therefore only the omp_level_master thread must work, i.e. : ! !$OMP BARRIER ! DO ind=1,ndomain ! IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE ! ... ! END DO ! !$OMP BARRIER ! There MUST be NO OMP BARRIER inside the DO-LOOP or any routine it calls. CONTAINS SUBROUTINE init_etat0 USE etat0_database_mod, ONLY: init_etat0_database => init_etat0 USE etat0_start_file_mod, ONLY: init_etat0_start_file => init_etat0 USE etat0_heldsz_mod, ONLY: init_etat0_held_suarez => init_etat0 CALL getin("etat0",etat0_type) SELECT CASE (TRIM(etat0_type)) CASE ('isothermal') CASE ('temperature_profile') CASE ('jablonowsky06') CASE ('dcmip5') CASE ('williamson91.6') CASE ('start_file') CALL init_etat0_start_file CASE ('database') CALL init_etat0_database CASE ('academic') CASE ('held_suarez') CALL init_etat0_held_suarez CASE ('venus') CASE ('dcmip1') CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') CASE ('dcmip3') CASE ('dcmip4') CASE ('dcmip2016_baroclinic_wave') CASE ('dcmip2016_cyclone') CASE ('dcmip2016_supercell') CASE ('bubble') CASE DEFAULT PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// & ' options are , , , , ,'& //' , , , , , ,' & //' , , ,'& //' , , ', 'bubble' STOP END SELECT END SUBROUTINE init_etat0 SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q) USE disvert_mod ! Generic interface USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0 USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0 USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0 USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0 USE etat0_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0 USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0 USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0 USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : getin_etat0_dcmip2016_baroclinic_wave=>getin_etat0 USE etat0_dcmip2016_cyclone_mod, ONLY : getin_etat0_dcmip2016_cyclone=>getin_etat0 USE etat0_dcmip2016_supercell_mod, ONLY : getin_etat0_dcmip2016_supercell=>getin_etat0 ! Ad hoc interfaces USE etat0_academic_mod, ONLY : etat0_academic=>etat0 USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0 USE etat0_venus_mod, ONLY : etat0_venus=>etat0 USE etat0_database_mod, ONLY : etat0_database=>etat0 USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0 TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_mass(:) TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_geopot(:) TYPE(t_field),POINTER :: f_w(:) TYPE(t_field),POINTER :: f_q(:) REAL(rstd),POINTER :: ps(:), mass(:,:) REAL :: etat0_ps_white_noise REAL :: etat0_theta_rhodz_white_noise REAL :: etat0_u_white_noise REAL :: etat0_q_white_noise LOGICAL :: autoinit_mass, collocated INTEGER :: ind ! most etat0 routines set ps and not mass ! in that case and if caldyn_eta == eta_lag ! the initial distribution of mass is taken to be the same ! as what the mass coordinate would dictate ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE. ! otherwise mass will be overwritten autoinit_mass = (caldyn_eta == eta_lag) etat0_type='jablonowsky06' CALL getin("etat0",etat0_type) !------------------- Generic interface --------------------- collocated=.TRUE. SELECT CASE (TRIM(etat0_type)) CASE ('isothermal') CALL getin_etat0_isothermal CASE ('temperature_profile') CALL getin_etat0_temperature CASE ('jablonowsky06') CASE ('dcmip1') CALL getin_etat0_dcmip1 CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') CALL getin_etat0_dcmip2 CASE ('dcmip3') CASE ('dcmip4') CALL getin_etat0_dcmip4 CASE ('dcmip5') CALL getin_etat0_dcmip5 CASE ('bubble') CALL getin_etat0_bubble CASE ('williamson91.6') autoinit_mass=.FALSE. CALL getin_etat0_williamson CASE ('dcmip2016_baroclinic_wave') CALL getin_etat0_dcmip2016_baroclinic_wave CASE ('dcmip2016_cyclone') CALL getin_etat0_dcmip2016_cyclone CASE ('dcmip2016_supercell') CALL getin_etat0_dcmip2016_supercell CASE DEFAULT collocated=.FALSE. autoinit_mass = .FALSE. END SELECT !------------------- Ad hoc interfaces -------------------- SELECT CASE (TRIM(etat0_type)) CASE ('database') CALL etat0_database(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('start_file') CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('academic') CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('held_suarez') PRINT *,"Held & Suarez (1994) test case" CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q) CASE ('venus') CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q) PRINT *, "Venus (Lebonnois et al., 2012) test case" CASE DEFAULT IF(collocated) THEN CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q) ELSE PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// & ' options are , , , , ,'& //' , , , , , ,' & //' , , ,'& //' , , ' STOP END IF END SELECT IF(autoinit_mass) THEN DO ind=1,ndomain IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) mass=f_mass(ind); ps=f_ps(ind) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps END DO END IF etat0_ps_white_noise=0. CALL getin("etat0_ps_white_noise",etat0_ps_white_noise) CALL add_white_noise(f_ps, etat0_ps_white_noise) etat0_theta_rhodz_white_noise=0. CALL getin("etat0_theta_rhodz_white_noise",etat0_theta_rhodz_white_noise) CALL add_white_noise(f_theta_rhodz, etat0_theta_rhodz_white_noise) etat0_u_white_noise=0. CALL getin("etat0_u_white_noise",etat0_u_white_noise) CALL add_white_noise(f_u, etat0_u_white_noise) etat0_q_white_noise=0. CALL getin("etat0_q_white_noise",etat0_q_white_noise) CALL add_white_noise(f_q, etat0_q_white_noise) END SUBROUTINE etat0 SUBROUTINE add_white_noise(field, factor) USE icosa IMPLICIT NONE TYPE(t_field),POINTER :: field(:) ! INOUT REAL,INTENT(IN) :: factor INTEGER,ALLOCATABLE :: seed(:) REAL,ALLOCATABLE :: random2d(:) REAL,ALLOCATABLE :: random3d(:,:) REAL,ALLOCATABLE :: random4d(:,:,:) REAL,POINTER :: field2d(:) REAL,POINTER :: field3d(:,:) REAL,POINTER :: field4d(:,:,:) INTEGER :: ind INTEGER :: m CALL RANDOM_SEED(SIZE=m) ALLOCATE(seed(m)) DO ind=1,ndomain IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE seed=domloc_glo_ind(ind) ! to be reproducible CALL RANDOM_SEED(PUT=seed) IF (field(ind)%ndim==2) THEN field2d=field(ind) ALLOCATE(random2d(size(field2d,1))) CALL RANDOM_NUMBER(random2d) field2d=field2d*(1.+random2d*factor) DEALLOCATE(random2d) ELSE IF (field(ind)%ndim==3) THEN field3d=field(ind) ALLOCATE(random3d(size(field3d,1),size(field3d,2))) CALL RANDOM_NUMBER(random3d) field3d=field3d*(1.+random3d*factor) DEALLOCATE(random3d) ELSE IF (field(ind)%ndim==4) THEN field4d=field(ind) ALLOCATE(random4d(size(field4d,1),size(field4d,2),size(field4d,3))) CALL RANDOM_NUMBER(random4d) field4d=field4d*(1.+random4d*factor) DEALLOCATE(random4d) ENDIF ENDDO END SUBROUTINE SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q) USE theta2theta_rhodz_mod TYPE(t_field),POINTER :: f_ps(:) TYPE(t_field),POINTER :: f_mass(:) TYPE(t_field),POINTER :: f_phis(:) TYPE(t_field),POINTER :: f_theta_rhodz(:) TYPE(t_field),POINTER :: f_u(:) TYPE(t_field),POINTER :: f_geopot(:) TYPE(t_field),POINTER :: f_W(:) TYPE(t_field),POINTER :: f_q(:) TYPE(t_field),POINTER,SAVE :: f_temp(:) REAL(rstd),POINTER :: ps(:) REAL(rstd),POINTER :: mass(:,:) REAL(rstd),POINTER :: phis(:) REAL(rstd),POINTER :: theta_rhodz(:,:,:) REAL(rstd),POINTER :: temp(:,:) REAL(rstd),POINTER :: u(:,:) REAL(rstd),POINTER :: geopot(:,:) REAL(rstd),POINTER :: W(:,:) REAL(rstd),POINTER :: q(:,:,:) INTEGER :: ind CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') DO ind=1,ndomain IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE ! IF (.NOT. assigned_domain(ind)) CYCLE CALL swap_dimensions(ind) CALL swap_geometry(ind) ps=f_ps(ind) mass=f_mass(ind) phis=f_phis(ind) theta_rhodz=f_theta_rhodz(ind) temp=f_temp(ind) u=f_u(ind) geopot=f_geopot(ind) w=f_w(ind) q=f_q(ind) IF( TRIM(etat0_type)=='williamson91.6' ) THEN CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q) ELSE CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q) ENDIF IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) ENDDO CALL deallocate_field(f_temp) END SUBROUTINE etat0_collocated SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) USE icosa USE pression_mod USE exner_mod USE omp_para REAL(rstd),INTENT(IN) :: ps(iim*jjm) REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm) REAL(rstd),INTENT(IN) :: q(iim*jjm,llm,nqtot) REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) INTEGER,INTENT(IN) :: offset REAL(rstd) :: p(iim*jjm,llm+1) REAL(rstd) :: cppd,Rd, mass, p_ij, chi,nu, entropy, theta INTEGER :: i,j,ij,l cppd=cpp Rd=kappa*cppd CALL compute_pression(ps,p,offset) ! flush p DO l = ll_begin, ll_end DO j=jj_begin-offset,jj_end+offset DO i=ii_begin-offset,ii_end+offset ij=(j-1)*iim+i mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass p_ij = .5*(p(ij,l)+p(ij,l+1)) ! pressure at full level SELECT CASE(caldyn_thermo) CASE(thermo_theta) theta = temp(ij,l)*(p_ij/preff)**(-kappa) theta_rhodz(ij,l) = mass * theta CASE(thermo_entropy) nu = log(p_ij/preff) chi = log(temp(ij,l)/Treff) entropy = cppd*chi-Rd*nu theta_rhodz(ij,l) = mass * entropy ! CASE(thermo_moist) ! q_ij=q(ij,l,1) ! r_ij=1.-q_ij ! mass=mass*(1-q_ij) ! dry mass ! nu = log(p_ij/preff) ! chi = log(temp(ij,l)/Treff) ! entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) ! theta_rhodz(ij,l) = mass * entropy CASE DEFAULT STOP END SELECT ENDDO ENDDO ENDDO END SUBROUTINE compute_temperature2entropy SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q) USE wind_mod USE disvert_mod USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0 USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0 USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0 USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0 USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0 USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0 USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0 USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0 USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0 REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm) REAL(rstd),INTENT(OUT) :: phis(iim*jjm) REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm) REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1) REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) REAL(rstd) :: ulon_i(iim*jjm,llm) REAL(rstd) :: ulat_i(iim*jjm,llm) REAL(rstd) :: ps_e(3*iim*jjm) REAL(rstd) :: mass_e(3*iim*jjm,llm) REAL(rstd) :: phis_e(3*iim*jjm) REAL(rstd) :: temp_e(3*iim*jjm,llm) REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) REAL(rstd) :: ulon_e(3*iim*jjm,llm) REAL(rstd) :: ulat_e(3*iim*jjm,llm) REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot) INTEGER :: l,ij REAL :: p_ik, v_ik, mass_ik LOGICAL :: autoinit_mass, autoinit_NH ! For NH geopotential and vertical momentum must be initialized. ! Unless autoinit_NH is set to .FALSE. , they will be initialized ! to hydrostatic geopotential and zero autoinit_mass = .TRUE. autoinit_NH = .NOT. hydrostatic w(:,:) = 0 SELECT CASE (TRIM(etat0_type)) CASE ('isothermal') CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE ('temperature_profile') CALL compute_etat0_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_etat0_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE('jablonowsky06') CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) CASE('dcmip1') CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) CASE('dcmip3') CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot CASE('dcmip4') CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip4(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE('dcmip5') CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE('bubble') CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) ! autoinit_NH = .FALSE. ! compute_bubble initializes geopot CASE('williamson91.6') CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1)) CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1)) autoinit_mass = .FALSE. ! do not overwrite mass CASE('dcmip2016_baroclinic_wave') CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE('dcmip2016_cyclone') CALL compute_dcmip2016_cyclone(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip2016_cyclone(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) CASE('dcmip2016_supercell') CALL compute_dcmip2016_supercell(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) CALL compute_dcmip2016_supercell(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) END SELECT IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps IF(autoinit_NH) THEN geopot(:,1) = phis(:) ! surface geopotential DO l = 1, llm DO ij=1,iim*jjm ! hybrid pressure coordinate p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g ! v=R.T/p, R=kappa*cpp v_ik = kappa*cpp*temp_i(ij,l)/p_ik geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g END DO END DO END IF CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u) END SUBROUTINE compute_etat0_collocated !----------------------------- Resting isothermal state -------------------------------- SUBROUTINE getin_etat0_isothermal etat0_temp=300 CALL getin("etat0_isothermal_temp",etat0_temp) END SUBROUTINE getin_etat0_isothermal SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q) INTEGER, INTENT(IN) :: ngrid REAL(rstd),INTENT(OUT) :: phis(ngrid) REAL(rstd),INTENT(OUT) :: ps(ngrid) REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) phis(:)=0 ps(:)=preff temp(:,:)=etat0_temp ulon(:,:)=0 ulat(:,:)=0 q(:,:,:)=0 END SUBROUTINE compute_etat0_isothermal END MODULE etat0_mod