Index: dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0.f90 =================================================================== --- dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0.f90 (revision 3861) +++ dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0.f90 (revision 3862) @@ -28,4 +28,5 @@ USE etat0_venus_mod, ONLY : etat0_venus=>etat0 USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0 + USE etat0_database_mod, ONLY : etat0_database=>etat0 IMPLICIT NONE @@ -72,4 +73,6 @@ CASE ('start_file') CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q) + CASE ('database') + CALL etat0_database(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) Index: dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0_database.f90 =================================================================== --- dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0_database.f90 (revision 3862) +++ dynamico_lmdz/aquaplanet/ICOSAGCM/src/etat0_database.f90 (revision 3862) @@ -0,0 +1,129 @@ +MODULE etat0_database_mod + + +CONTAINS + + SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q) + USE icosa + USE restart_mod + USE wind_mod + USE write_field_mod + USE time_mod + USE transfert_mod + USE xios_mod + USE write_field_mod + USE vertical_remap_mod + USE theta2theta_rhodz_mod + USE qsat_mod + USE pression_mod + USE write_etat0_mod + IMPLICIT NONE + TYPE(t_field),POINTER :: f_ps(:) + TYPE(t_field),POINTER :: f_phis(:) + TYPE(t_field),POINTER :: f_theta_rhodz(:) + TYPE(t_field),POINTER :: f_u(:) + TYPE(t_field),POINTER :: f_q(:) + + TYPE(t_field),POINTER :: f_ulon_reg(:) + TYPE(t_field),POINTER :: f_ulat_reg(:) + TYPE(t_field),POINTER :: f_temp_reg(:) + TYPE(t_field),POINTER :: f_q_reg(:) + + TYPE(t_field),POINTER :: f_ts(:) + TYPE(t_field),POINTER :: f_z(:) + TYPE(t_field),POINTER :: f_ulon(:) + TYPE(t_field),POINTER :: f_ulat(:) + TYPE(t_field),POINTER :: f_temp(:) + TYPE(t_field),POINTER :: f_q1(:) + TYPE(t_field),POINTER :: f_qsat(:) + TYPE(t_field),POINTER :: f_p(:) + INTEGER :: nb_level + REAL,ALLOCATABLE:: levels(:) + INTEGER :: ind + + CALL xios_read_field("relief",f_phis) + CALL writeField("relief_out",f_phis,once=.TRUE.) + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind)) CYCLE + f_phis(ind)%rval2d(:)=f_phis(ind)%rval2d(:)*g + ENDDO + + + CALL xios_get_axis_attr("lev_ecdyn",n_glo=nb_level) + ALLOCATE(levels(nb_level)) + CALL xios_get_axis_attr("lev_ecdyn",value=levels) + levels=levels*100 ! hectoPascal -> Pascal + + CALL allocate_field(f_ts, field_t, type_real, name="ts") + CALL allocate_field(f_z, field_t, type_real, name="z") + CALL allocate_field(f_ulon_reg, field_t, type_real,nb_level) + CALL allocate_field(f_ulat_reg, field_t, type_real,nb_level) + CALL allocate_field(f_temp_reg, field_t, type_real,nb_level) + CALL allocate_field(f_q_reg, field_t, type_real,nb_level) + + CALL allocate_field(f_q1, field_t, type_real,llm) + CALL allocate_field(f_qsat, field_t, type_real,llm) + CALL allocate_field(f_p, field_t, type_real,llm+1) + CALL allocate_field(f_temp, field_t, type_real,llm) + CALL allocate_field(f_ulon, field_t, type_real,llm) + CALL allocate_field(f_ulat, field_t, type_real,llm) + + CALL xios_read_field("z",f_z) + CALL xios_read_field("ps",f_ps) + CALL xios_read_field("ts",f_ts) + CALL writeField("ps_out",f_ps) +! CALL writeField("phis_out",f_phis,once=.TRUE.) +! CALL writeField("ts_out",f_ts,once=.TRUE.) + +! make correction to ps due to relief at higher resolution +! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind)) CYCLE + f_ps(ind)%rval2d(:)=f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:)-f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:)) + ENDDO + CALL transfert_request(f_ps,req_i0) + CALL writeField("ps_out",f_ps) + + + + CALL xios_read_field("temp",f_temp_reg) + CALL vertical_remap(levels,f_temp_reg,f_ps,f_temp) + CALL transfert_request(f_temp,req_i0) + CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz) + + CALL xios_read_field("u",f_ulon_reg) + CALL vertical_remap(levels,f_ulon_reg,f_ps,f_ulon) + CALL xios_read_field("v",f_ulat_reg) + CALL vertical_remap(levels,f_ulat_reg,f_ps,f_ulat) + CALL transfert_request(f_ulat,req_i0) + CALL transfert_request(f_ulon,req_i0) + CALL ulonlat2un(f_ulon, f_ulat,f_u) + + CALL xios_read_field("q",f_q_reg) + CALL vertical_remap(levels,f_q_reg,f_ps,f_q1) + + CALL pression(f_ps,f_p) +! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2 + CALL qsat(f_temp,f_p,f_qsat) + CALL transfert_request(f_qsat,req_i0) + + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind)) CYCLE + f_q(ind)%rval4d(:,:,:)=0 + f_q(ind)%rval4d(:,:,1)=f_q1(ind)%rval3d(:,:)*f_qsat(ind)%rval3d(:,:)*0.01 + WHERE(f_q(ind)%rval4d(:,:,1)<0) f_q(ind)%rval4d(:,:,1)=0 + ENDDO + + + CALL writeField("ulon_out",f_ulon,once=.TRUE.) + CALL writeField("ulat_out",f_ulat,once=.TRUE.) + CALL writeField("temp_out",f_temp,once=.TRUE.) + CALL writeField("temp2_out",f_theta_rhodz,once=.TRUE.) + CALL writeField("f_qsat",f_qsat,once=.TRUE.) + CALL writeField("f_q1",f_q1,once=.TRUE.) + CALL writeField("f_q",f_q,once=.TRUE.) + CALL write_etat0(0,f_ps,f_phis,f_theta_rhodz,f_u, f_q) + + END SUBROUTINE etat0 + +END MODULE etat0_database_mod Index: dynamico_lmdz/aquaplanet/ICOSAGCM/src/q_sat.f90 =================================================================== --- dynamico_lmdz/aquaplanet/ICOSAGCM/src/q_sat.f90 (revision 3862) +++ dynamico_lmdz/aquaplanet/ICOSAGCM/src/q_sat.f90 (revision 3862) @@ -0,0 +1,98 @@ +MODULE qsat_mod + + + + +CONTAINS + + SUBROUTINE qsat(f_temp,f_p,f_qs) + USE icosa + IMPLICIT NONE + TYPE(t_field), POINTER :: f_temp(:) ! IN : temperature + TYPE(t_field), POINTER :: f_p(:) ! IN : pressure at mid-levels + TYPE(t_field), POINTER :: f_qs(:) ! OUT : vapeur d'eau saturante en kg/kg + + REAL(rstd),POINTER :: temp(:,:), p(:,:), qs(:,:) + INTEGER :: ind + + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind)) CYCLE + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + temp=f_temp(ind) + p=f_p(ind) + qs=f_qs(ind) + CALL compute_qsat(temp,p,qs) + END DO + + END SUBROUTINE qsat + + + SUBROUTINE compute_qsat(temp,p,qsat) + USE icosa + USE omp_para + IMPLICIT NONE + +!====================================================================== +! Autheur(s): Z.X. Li (LMD/CNRS) +! reecriture vectorisee par F. Hourdin. +! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.) +!====================================================================== +! Arguments: +! kelvin---input-R: temperature en Kelvin +! millibar--input-R: pression en mb +! +! q_sat----output-R: vapeur d'eau saturante en kg/kg +!====================================================================== +! + REAL,INTENT(IN) :: temp(iim*jjm,llm) + REAL,INTENT(IN) :: p (iim*jjm,llm+1) + REAL,INTENT(OUT) :: qsat(iim*jjm,llm) + + REAL, PARAMETER :: r2es=611.14 *18.0153/28.9644 + + REAL :: r3es + REAL, PARAMETER :: r3les=17.269 + REAL, PARAMETER :: r3ies=21.875 + + REAL :: r4es + REAL, PARAMETER :: r4les=35.86 + REAL, PARAMETER :: r4ies=7.66 + + REAL, PARAMETER :: rtt=273.16 + REAL, PARAMETER :: retv=28.9644/18.0153 - 1.0 + + REAL :: zqsat, pmid + INTEGER :: l,i,j,ij +! +! ------------------------------------------------------------------ +! +! + + DO l=ll_begin,ll_end + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + + IF (temp(ij,l) .LE. rtt) THEN + r3es = r3ies + r4es = r4ies + ELSE + r3es = r3les + r4es = r4les + ENDIF + pmid=0.5*(p(ij,l)+p(ij,l+1)) + zqsat=r2es/pmid*EXP(r3es*(temp(ij,l)-rtt)/(temp(ij,l)-r4es)) + zqsat=MIN(0.5,zqsat) + zqsat=zqsat/(1.-retv *zqsat) + + qsat(ij,l)= zqsat + + ENDDO + ENDDO + ENDDO + + + END SUBROUTINE compute_qsat + +END MODULE qsat_mod Index: dynamico_lmdz/aquaplanet/ICOSAGCM/src/vertical_remap.f90 =================================================================== --- dynamico_lmdz/aquaplanet/ICOSAGCM/src/vertical_remap.f90 (revision 3862) +++ dynamico_lmdz/aquaplanet/ICOSAGCM/src/vertical_remap.f90 (revision 3862) @@ -0,0 +1,110 @@ +MODULE vertical_remap_mod + USE icosa + + +CONTAINS + + + SUBROUTINE vertical_remap(pressure_level,field_in,f_ps,field_out) + USE icosa + USE pression_mod + USE omp_para + IMPLICIT NONE + REAL(rstd), INTENT(IN) :: pressure_level(:) + TYPE(t_field),POINTER :: field_in(:) + TYPE(t_field),POINTER :: f_ps(:) + TYPE(t_field),POINTER :: field_out(:) + + TYPE(t_field),POINTER :: f_p(:) + REAL(rstd),POINTER :: in(:,:) + REAL(rstd),POINTER :: out(:,:) + REAL(rstd),POINTER :: p(:,:) + + INTEGER :: ind + + CALL allocate_field(f_p,field_t,type_real,llm+1) + CALL pression(f_ps,f_p) + + DO ind=1,ndomain + IF (.NOT. assigned_domain(ind)) CYCLE + CALL swap_dimensions(ind) + CALL swap_geometry(ind) + p=f_p(ind) + in=field_in(ind) + out=field_out(ind) + CALL compute_vertical_remap(pressure_level,in,p,out) + ENDDO + + END SUBROUTINE vertical_remap + + SUBROUTINE compute_vertical_remap(pressure_level,in,p,out) + USE omp_para + IMPLICIT NONE + REAL(rstd),INTENT(IN) :: pressure_level(:) + REAL(rstd),INTENT(IN) :: in(:,:) + REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1) + REAL(rstd),INTENT(OUT) :: out(iim*jjm,llm) + REAL(rstd) :: coeff, pmid + INTEGER :: i,j,ij,l,n,nb_level + INTEGER :: a + INTEGER :: b + LOGICAL :: positive + + nb_level=size(pressure_level) + IF (pressure_level(1)>=pressure_level(nb_level)) THEN + positive=.FALSE. + ELSE + positive=.TRUE. + ENDIF + + !$OMP BARRIER + IF (is_omp_level_master) THEN + + DO l=1,llm + DO j=jj_begin,jj_end + DO i=ii_begin,ii_end + ij=(j-1)*iim+i + pmid=0.5*(p(ij,l)+p(ij,l+1)) + IF (positive) THEN + a=0 + DO n=1,nb_level-1 + IF ( (pmid>=pressure_level(n) .AND. pmidpressure_level(n+1))) THEN + a=n ; b=n+1 ; EXIT + ENDIF + ENDDO + + IF (a==0) THEN + IF (pmid >= pressure_level(1)) THEN + a=1 ; b=2 + ELSE + a=nb_level-1 ; b=nb_level + ENDIF + ENDIF + ENDIF + + coeff=(pmid-pressure_level(a))/(pressure_level(a)-pressure_level(b)) + out(ij,l)=in(ij,a)+coeff*(in(ij,a)-in(ij,b)) + ENDDO + ENDDO + ENDDO + + ENDIF + !$OMP BARRIER + + END SUBROUTINE compute_vertical_remap + +END MODULE vertical_remap_mod