Changeset 847


Ignore:
Timestamp:
05/05/19 00:08:28 (5 years ago)
Author:
dubos
Message:

devel : split etat0 into etat0, etat0_collocated and etat0_isothermal

Location:
codes/icosagcm/devel/src/initial
Files:
2 added
1 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/devel/src/initial/etat0.f90

    r732 r847  
    11MODULE etat0_mod 
    22  USE icosa 
    3   USE omp_para 
    4   USE caldyn_vars_mod 
     3  USE etat0_collocated_mod 
    54  IMPLICIT NONE          
    65  PRIVATE 
    76 
    8     CHARACTER(len=255),SAVE :: etat0_type 
    9 !$OMP THREADPRIVATE(etat0_type) 
    10  
    11     REAL(rstd) :: etat0_temp 
    12  
    13     PUBLIC :: etat0, init_etat0, etat0_type 
     7  PUBLIC :: etat0, init_etat0 
    148 
    159! Important notes for OpenMP 
     
    6761 
    6862  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q) 
     63    USE omp_para, ONLY : is_omp_level_master 
    6964    USE mpipara, ONLY : is_mpi_root 
    7065    USE disvert_mod 
    7166    ! Generic interface 
     67    USE etat0_isothermal_mod, ONLY : getin_etat0_isothermal=>getin_etat0 
    7268    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0 
    7369    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0 
     
    183179  END SUBROUTINE etat0 
    184180 
    185   SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q) 
    186     USE theta2theta_rhodz_mod 
    187     TYPE(t_field),POINTER :: f_ps(:) 
    188     TYPE(t_field),POINTER :: f_mass(:) 
    189     TYPE(t_field),POINTER :: f_phis(:) 
    190     TYPE(t_field),POINTER :: f_theta_rhodz(:) 
    191     TYPE(t_field),POINTER :: f_u(:) 
    192     TYPE(t_field),POINTER :: f_geopot(:) 
    193     TYPE(t_field),POINTER :: f_W(:) 
    194     TYPE(t_field),POINTER :: f_q(:) 
    195    
    196     TYPE(t_field),POINTER,SAVE :: f_temp(:) 
    197     REAL(rstd),POINTER :: ps(:) 
    198     REAL(rstd),POINTER :: mass(:,:) 
    199     REAL(rstd),POINTER :: phis(:) 
    200     REAL(rstd),POINTER :: theta_rhodz(:,:,:) 
    201     REAL(rstd),POINTER :: temp(:,:) 
    202     REAL(rstd),POINTER :: u(:,:) 
    203     REAL(rstd),POINTER :: geopot(:,:) 
    204     REAL(rstd),POINTER :: W(:,:) 
    205     REAL(rstd),POINTER :: q(:,:,:) 
    206     INTEGER :: ind 
    207  
    208     CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 
    209  
    210     DO ind=1,ndomain 
    211       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE 
    212 !      IF (.NOT. assigned_domain(ind)) CYCLE 
    213       CALL swap_dimensions(ind) 
    214       CALL swap_geometry(ind) 
    215       ps=f_ps(ind) 
    216       mass=f_mass(ind) 
    217       phis=f_phis(ind) 
    218       theta_rhodz=f_theta_rhodz(ind) 
    219       temp=f_temp(ind) 
    220       u=f_u(ind) 
    221       geopot=f_geopot(ind) 
    222       w=f_w(ind) 
    223       q=f_q(ind) 
    224  
    225       IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
    226          CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q) 
    227       ELSE 
    228          CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q) 
    229       ENDIF 
    230  
    231       IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1) 
    232      
    233     ENDDO 
    234      
    235     CALL deallocate_field(f_temp) 
    236      
    237   END SUBROUTINE etat0_collocated 
    238  
    239   SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset) 
    240     USE icosa 
    241     USE pression_mod 
    242     USE exner_mod 
    243     USE omp_para 
    244     REAL(rstd),INTENT(IN)  :: ps(iim*jjm) 
    245     REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm) 
    246     REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot) 
    247     REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm) 
    248     INTEGER,INTENT(IN) :: offset 
    249  
    250     REAL(rstd) :: p(iim*jjm,llm+1) 
    251     REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta 
    252     INTEGER :: i,j,ij,l 
    253  
    254     cppd=cpp 
    255     Rd=kappa*cppd 
    256  
    257     CALL compute_pression(ps,p,offset) 
    258     ! flush p 
    259     DO    l    = ll_begin, ll_end 
    260        DO j=jj_begin-offset,jj_end+offset 
    261           DO i=ii_begin-offset,ii_end+offset 
    262              ij=(j-1)*iim+i 
    263              mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass 
    264              p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level 
    265              SELECT CASE(caldyn_thermo) 
    266              CASE(thermo_theta) 
    267                 theta = temp(ij,l)*(p_ij/preff)**(-kappa)  
    268                 theta_rhodz(ij,l) = mass * theta 
    269              CASE(thermo_entropy) 
    270                 nu = log(p_ij/preff) 
    271                 chi = log(temp(ij,l)/Treff) 
    272                 entropy = cppd*chi-Rd*nu 
    273                 theta_rhodz(ij,l) = mass * entropy 
    274 !             CASE(thermo_moist) 
    275 !                q_ij=q(ij,l,1) 
    276 !                r_ij=1.-q_ij 
    277 !                mass=mass*(1-q_ij) ! dry mass 
    278 !                nu = log(p_ij/preff) 
    279 !                chi = log(temp(ij,l)/Treff) 
    280 !                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu) 
    281 !                theta_rhodz(ij,l) = mass * entropy                 
    282                 CASE DEFAULT 
    283                    STOP 
    284              END SELECT 
    285           ENDDO 
    286        ENDDO 
    287     ENDDO 
    288   END SUBROUTINE compute_temperature2entropy 
    289  
    290   SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q) 
    291     USE wind_mod 
    292     USE disvert_mod 
    293     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0 
    294     USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0 
    295     USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0 
    296     USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0 
    297     USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0 
    298     USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0 
    299     USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0   
    300     USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0 
    301     USE etat0_temperature_mod, ONLY: compute_etat0_temperature => compute_etat0 
    302     USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0 
    303     USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0 
    304     USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0 
    305     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm) 
    306     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm) 
    307     REAL(rstd),INTENT(OUT) :: phis(iim*jjm) 
    308     REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm) 
    309     REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm) 
    310     REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1) 
    311     REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1) 
    312     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot) 
    313  
    314     REAL(rstd) :: ulon_i(iim*jjm,llm) 
    315     REAL(rstd) :: ulat_i(iim*jjm,llm) 
    316  
    317     REAL(rstd) :: ps_e(3*iim*jjm) 
    318     REAL(rstd) :: mass_e(3*iim*jjm,llm) 
    319     REAL(rstd) :: phis_e(3*iim*jjm) 
    320     REAL(rstd) :: temp_e(3*iim*jjm,llm) 
    321     REAL(rstd) :: geopot_e(3*iim*jjm,llm+1) 
    322     REAL(rstd) :: ulon_e(3*iim*jjm,llm) 
    323     REAL(rstd) :: ulat_e(3*iim*jjm,llm) 
    324     REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot) 
    325  
    326     INTEGER :: l,i,j,ij 
    327     REAL :: p_ik, v_ik, mass_ik 
    328     LOGICAL :: autoinit_mass, autoinit_NH 
    329  
    330     ! For NH geopotential and vertical momentum must be initialized. 
    331     ! Unless autoinit_NH is set to .FALSE. , they will be initialized 
    332     ! to hydrostatic geopotential and zero 
    333     autoinit_mass = .TRUE. 
    334     autoinit_NH = .NOT. hydrostatic 
    335     w(:,:) = 0 
    336  
    337     SELECT CASE (TRIM(etat0_type)) 
    338     CASE ('isothermal') 
    339        CALL compute_etat0_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q) 
    340        CALL compute_etat0_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    341     CASE ('temperature_profile') 
    342        CALL compute_etat0_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q) 
    343        CALL compute_etat0_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    344     CASE('jablonowsky06') 
    345        CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 
    346        CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e) 
    347     CASE('dcmip1') 
    348        CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    349        CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    350     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear') 
    351        CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i) 
    352        CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)       
    353     CASE('dcmip3') 
    354        CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) 
    355        CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 
    356        autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot 
    357     CASE('dcmip4') 
    358        CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    359        CALL compute_dcmip4(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    360     CASE('dcmip5') 
    361        CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    362        CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    363     CASE('bubble') 
    364        CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q) 
    365        CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e) 
    366 !       autoinit_NH = .FALSE. ! compute_bubble initializes geopot 
    367     CASE('williamson91.6') 
    368        CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1)) 
    369        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)) 
    370        autoinit_mass = .FALSE. ! do not overwrite mass 
    371     CASE('dcmip2016_baroclinic_wave') 
    372        CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    373        CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    374     CASE('dcmip2016_cyclone') 
    375        CALL compute_dcmip2016_cyclone(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    376        CALL compute_dcmip2016_cyclone(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    377     CASE('dcmip2016_supercell') 
    378        CALL compute_dcmip2016_supercell(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q) 
    379        CALL compute_dcmip2016_supercell(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e) 
    380     END SELECT 
    381  
    382     IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps 
    383     IF(autoinit_NH) THEN 
    384        geopot(:,1) = phis(:) ! surface geopotential 
    385        DO l = 1, llm 
    386           DO ij=1,iim*jjm 
    387              ! hybrid pressure coordinate 
    388              p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) 
    389              mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g 
    390              ! v=R.T/p, R=kappa*cpp 
    391              v_ik = kappa*cpp*temp_i(ij,l)/p_ik 
    392              geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g 
    393           END DO 
    394        END DO 
    395     END IF 
    396  
    397     CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u) 
    398  
    399   END SUBROUTINE compute_etat0_collocated 
    400  
    401 !----------------------------- Resting isothermal state -------------------------------- 
    402  
    403   SUBROUTINE getin_etat0_isothermal 
    404     etat0_temp=300 
    405     CALL getin("etat0_isothermal_temp",etat0_temp) 
    406   END SUBROUTINE getin_etat0_isothermal 
    407  
    408   SUBROUTINE compute_etat0_isothermal(ngrid, phis, ps, temp, ulon, ulat, q) 
    409     INTEGER, INTENT(IN)    :: ngrid 
    410     REAL(rstd),INTENT(OUT) :: phis(ngrid) 
    411     REAL(rstd),INTENT(OUT) :: ps(ngrid) 
    412     REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) 
    413     REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) 
    414     REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) 
    415     REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) 
    416     phis(:)=0 
    417     ps(:)=preff 
    418     temp(:,:)=etat0_temp 
    419     ulon(:,:)=0 
    420     ulat(:,:)=0 
    421     q(:,:,:)=0 
    422   END SUBROUTINE compute_etat0_isothermal 
    423  
    424181END MODULE etat0_mod 
Note: See TracChangeset for help on using the changeset viewer.