MODULE wrf_lmdz_mod ! ! Module for the interface to include LMDZ (http://lmdz.lmd.jussieu.fr/) physical packages in WRF ! ! L. Fita, Laboratoire Meterologie Dynamique, IPSL, UPMC, CNRS, Jussieu, Paris, France ! November 2013 ! ! WRF modules ! USE module_domain_type !!!!!!! Subroutines ! wrf_varKsoil: Subroutine to provide WRF variable values from the 4 LMDZ soil types ! lmdz_varKsoil: Subroutine to provide variable values to the 4 LMDZ soil types according to different methods ! LMDZmoist_WRFmoist: Subroutine to transform from LMDZ moisture array to the WRF moisture array ! WRFmoist_LMDZmoist: Subroutine to transform from WRF moisture array to the LMDZ moisture array ! eta_to_pressure: Subroutine to transform from eta levels to pressure levels ! temp_to_theta: Subroutine to transform from temperatures to WRF potential temperatures ! theta_to_temp: Subroutine to transform from WRF potential temperatures to temperature ! temp_to_theta1: Subroutine to transform from temperature to potential temperature ! theta_to_temp1: Subroutine to transform from potential temperature to temperature ! def_get_scalar_value: Subroutine to get a scalar value from a .def file ! NOTE: For that scalar values with scientific notation (EN, or e-N) it will not work... ! next version! ! value_kind: Subroutine to determine which kind of value is given ! string_split: Subtoutine to split a string according to a character and returns the 'Nsec' of the split ! load_lmdz: Subroutine to load LMDZ variables with values from WRF ! load_lmdz_limit: Subroutine to load LMDZ limit variables with values from WRF ! get_lmdz: Subroutine to get LMDZ values and fill WRF variables ! get_lowbdy: Subroutine to load LMDZ limit variables with values from lowbdy WRF ! vect_mat: Subroutine to transform a LMDZ physic 1D vector to a 3D matrix ! vect_mat_zstagg: Subroutine to transform a LMDZ physic 1D vector to a 3D matrix ! mat_vect: Subroutine to transform a 3D matrix to a LMDZ physic 1D vector ! mat_vect_zstagg: Subroutine to transform a 3D matrix z-staggered to a LMDZ physic 1D vector ! WRFlanduse_LMDZKsoil: Subroutine to retrieve LMDZ Ksoil classification using WRF landuse ! interpolate_layers: Subroutine to interpolate values between different layers ! (constant values along the depth of a layer) ! interpolate_lin: Program to interpolate values y=a+b*x with: (from clWRF modifications) ! table_characteristics: Subroutine to read values from a WRF .TBL ! read_table: Subroutine to read values from a WRF .TBL ! compute_soil_mass: Subroutine to compute the total soil mass (kg) from a point with ! different percentages of soil types ! compute_TOTsoil_mass_water: Subroutine to compute total soil mass and water from a ! given grid point using a given data-set from SOILPARM.TBL ! compute_TOTsoil_mass_water_values: Subroutine to compute total soil mass and water ! from a given grid point using the values from a given data-set of SOILPARM.TBL ! get_lmdz_out: Subroutine to get all the LMDZ output variables into the WRF ! Registry structure ! put_lmdz_in_WRFout: Subroutine to get LMDZ output and get into WRF standard output variables ! LMDZ modules ! USE indice_sol_mod ! IMPLICIT NONE ! INCLUDE "dimsoil.h" CONTAINS SUBROUTINE wrf_varKsoil(is,ie,js,je,dx,dy,wbdy,kmethod,ter,lic,oce,sic,vter,vlic, & voce,vsic,var) ! Subroutine to provide WRF variable values from the 4 LMDZ soil types IMPLICIT NONE INTEGER, INTENT(IN) :: is,ie,js,je,dx,dy,wbdy CHARACTER(LEN=50), INTENT(IN) :: kmethod REAL, DIMENSION(is:ie,js:je), INTENT(IN) :: ter,lic,oce,sic REAL, DIMENSION(is:ie,js:je), INTENT(IN) :: vter,vlic,voce,vsic REAL, DIMENSION(is:ie,js:je), INTENT(OUT) :: var ! Local INTEGER :: i,j,ddx,ddy CHARACTER(LEN=50) :: errmsg, fname REAL*8, DIMENSION(is:ie,js:je) :: prop !!!!!!! Variables ! d[x/y]: dimension of the matrices ! kmethod: method to use ! 'prod': product, variable is distributed by percentage of the soil type vsoil=var*soil ! ter: percentage of continent ! lic: percentage of frozen continent ! oce: percentage of ocean ! sic: percentage of frozen ocean ! vter: variable on continent ! vlic: variable on frozen continent ! voce: variable on ocean ! vsic: variable on frozen ocean ! prop: proportion between soil fractions ! var: variable to be reinterpreted from LMDZ soil types fname="'wrf_varksoil'" errmsg='ERROR -- error -- ERROR -- error' prop=-1.d0 var=0. ddx=dx-2*wbdy ddy=dy-2*wbdy methodkind: IF (TRIM(kmethod) == 'prod') THEN DO i=1,ddx DO j=1,ddy var(i,j)=vter(i,j)*ter(i,j)+vlic(i,j)*lic(i,j)+voce(i,j)*oce(i,j)+vsic(i,j)*sic(i,j) END DO END DO ELSE PRINT *,TRIM(errmsg) PRINT *,' ' // TRIM(fname) // ": method '" // TRIM(kmethod) // "' does not exist !!!!!" STOP END IF methodkind END SUBROUTINE wrf_varKsoil SUBROUTINE lmdz_varKsoil(is,ie,js,je,dx,dy,wbdy,var,kmethod,ter,lic,oce,sic,vter, & vlic,voce,vsic) ! Subroutine to provide variable values to the 4 LMDZ soil types according to ! different methods IMPLICIT NONE INTEGER, INTENT(IN) :: is,ie,js,je,dx,dy,wbdy REAL, DIMENSION(is:ie,js:je), INTENT(IN) :: var CHARACTER(LEN=50), INTENT(IN) :: kmethod REAL, DIMENSION(is:ie,js:je), INTENT(IN) :: ter,lic,oce,sic REAL, DIMENSION(is:ie,js:je), INTENT(OUT) :: vter,vlic,voce,vsic ! Local INTEGER :: i,j,ddx,ddy CHARACTER(LEN=50) :: errmsg, fname REAL*8, DIMENSION(is:ie,js:je) :: prop !!!!!!! Variables ! d[x/y]: dimension of the matrices ! var: variable to be reinterpreted in LMDZ soil types ! kmethod: method to use ! 'NOlic': NO lic, variable is not located in ice continent soil type: lic ! 'NOlnd': NO land, variable is distributed only to sea soil types: oce, sic ! 'NOneg': NO negative, variable is distributed only to free ice soil types: ter, oce ! 'NOoce': NO oce, variable is not located in ocean soil type: oce ! 'NOpos': NO positive, variable is distributed only to iced soil types: lic, sic ! 'NOsea': NO sea, variable is distributed only to land soil types: ter, lic ! 'NOsic': NO sic, variable is not located in ice ocean soil type: sic ! 'NOter': NO ter, variable is not located in land soil type: ter ! 'prod': product, variable is distributed by percentage of the soil type vsoil=var*soil ! 'lic:' lic, variable only is located in ice land soil type: lic ! 'oce:' oce, variable only is located in ocean soil type: oce ! 'sic:' sic, variable only is located in ice ocean soil type: sic ! 'ter:' ter, variable only is located in land soil type: ter ! ter: percentage of continent ! lic: percentage of frozen continent ! oce: percentage of ocean ! sic: percentage of frozen ocean ! vter: variable on continent ! vlic: variable on frozen continent ! voce: variable on ocean ! vsic: variable on frozen ocean ! prop: proportion between soil fractions fname="'lmdz_varksoil'" errmsg='ERROR -- error -- ERROR -- error' ddx=dx-2*wbdy ddy=dy-2*wbdy prop=-1.d0 vter=0. vlic=0. voce=0. vsic=0. methodkind: IF (TRIM(kmethod) == 'NOlic') THEN WHERE (ter+oce+sic /= 0.) prop=1.d0/(ter+oce+sic) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vter(i,j)=var(i,j)*ter(i,j)*prop(i,j) voce(i,j)=var(i,j)*oce(i,j)*prop(i,j) vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOlnd') THEN WHERE (oce+sic /= 0.) prop=oce*1.d0/(oce+sic) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN voce(i,j)=var(i,j)*prop(i,j) vsic(i,j)=var(i,j)*(1.-prop(i,j)) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOneg') THEN WHERE (ter+oce /= 0.) prop=ter*1.d0/(ter+oce) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vter(i,j)=var(i,j)*prop(i,j) voce(i,j)=var(i,j)*(1.-prop(i,j)) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOoce') THEN WHERE (ter+lic+sic /= 0.) prop=1.d0/(ter+lic+sic) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vter(i,j)=var(i,j)*ter(i,j)*prop(i,j) vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j) vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOpos') THEN WHERE (lic+sic /=0.) prop=lic*1.d0/(lic+sic) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vlic(i,j)=var(i,j)*prop(i,j) vsic(i,j)=var(i,j)*(1.-prop(i,j)) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOsea') THEN WHERE (ter+lic /= 0.) prop=ter*1.d0/(ter+lic) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vter(i,j)=var(i,j)*prop(i,j) vlic(i,j)=var(i,j)*(1.-prop(i,j)) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOsic') THEN WHERE (ter+lic+oce /= 0.) prop=1.d0/(ter+lic+oce) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vter(i,j)=var(i,j)*ter(i,j)*prop(i,j) vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j) voce(i,j)=var(i,j)*oce(i,j)*prop(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'NOter') THEN WHERE (lic+sic+oce /= 0.) prop=1.d0/(lic+sic+oce) DO i=1,ddx DO j=1,ddy IF (prop(i,j) >= 0.) THEN vlic(i,j)=var(i,j)*lic(i,j)*prop(i,j) vsic(i,j)=var(i,j)*sic(i,j)*prop(i,j) voce(i,j)=var(i,j)*oce(i,j)*prop(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'lic') THEN DO i=1,ddx DO j=1,ddy IF (lic(i,j) /= 0.) THEN vlic(i,j)=var(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'oce') THEN DO i=1,ddx DO j=1,ddy IF (oce(i,j) /= 0.) THEN voce(i,j)=var(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'prod') THEN DO i=1,ddx DO j=1,ddy vter(i,j)=var(i,j)*ter(i,j)*1.d0 vlic(i,j)=var(i,j)*lic(i,j)*1.d0 voce(i,j)=var(i,j)*oce(i,j)*1.d0 vsic(i,j)=var(i,j)*sic(i,j)*1.d0 END DO END DO ELSEIF (TRIM(kmethod) == 'sic') THEN DO i=1,ddx DO j=1,ddy IF (sic(i,j) /= 0.) THEN vsic(i,j)=var(i,j) END IF END DO END DO ELSEIF (TRIM(kmethod) == 'ter') THEN DO i=1,ddx DO j=1,ddy IF (ter(i,j) /= 0.) THEN vter(i,j)=var(i,j) END IF END DO END DO ELSE PRINT *,TRIM(errmsg) PRINT *,' ' // TRIM(fname) // ": method '" // TRIM(kmethod) // "' does not exist !!!!!" STOP END IF methodkind END SUBROUTINE lmdz_varKsoil SUBROUTINE LMDZmoist_WRFmoist(lmdzMixRat, dx, dy, dz, Nwmoist, wbdy, wqvid, wqcid, & & wqrid, Nlmdzq, wMOIST) ! Subroutine to transform from LMDZ moisture array to the WRF moisture array IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy INTEGER, INTENT(IN) :: Nwmoist, wqvid, wqcid, & wqrid INTEGER, INTENT(IN) :: Nlmdzq REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz,Nlmdzq), INTENT(IN) :: & lmdzMixRat REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1,Nwmoist), INTENT(OUT) :: wMOIST ! Local INTEGER :: i,j,k, iq INTEGER :: ddx, ddy INTEGER, ALLOCATABLE, DIMENSION(:) :: widlmdz REAL, DIMENSION(1-wbdy:dx-wbdy+1,1-wbdy:dy-wbdy+1) :: qcr_ratio, wcrMOIST !!!!!!! Variables ! widlmdz: vector with the id of the WRF tracers to be also included in the LMDZ moist array ! d[x,y,z]: Spatial dimension of the WRF moisture array ! Nwmoist: number of species in wMOIST ! wqvid: index for the water vapor mixing ratio ! wqcid: index for the cloud water mixing ratio ! wqrid: index for the rain water mixing ratio ! Nlmdzq: number of species in LMDZ ! lmdxMixRat: LMDZ mixing ratio array ! qcr_ratio: ratio between cloud/rain liquid water in WRF moisture array ! wcrMOIST: WRF array combination of cloud and rain liquid water ! wMOIST: WRF moisture array !!!!!! L. Fita, July 203 ! At this stage, we only care about water vapor, cloud liquid water and rain liquid water. ddx=dx-2*wbdy ddy=dy-2*wbdy wMOIST=0. wcrMOIST=0. DO k=1,dz DO j=1,ddy DO i=1,ddx wMOIST(i,k,j,wqvid) = lmdzMixRat(ddx*(j-1)+i,k,1) wMOIST(i,k,j,wqcid) = lmdzMixRat(ddx*(j-1)+i,k,2) !! wcrMOIST(i,j) = lmdzMixRat(ddx*(j-1)+1,k,2) END DO END DO !! L. Fita, July 2013 !! We can not know how to split liquid water between cloud and rain. Thus, we use the previous !! time-step proportion !! qcr_ratio=wMOIST(:,k,:,wqcid)/(wMOIST(:,k,:,wqcid) + wMOIST(:,k,:,wqrid)) !! wMOIST(:,k,:,wqcid) = wcrMOIST*qcr_ratio !! wMOIST(:,k,:,wqrid) = wcrMOIST*(1.-qcr_ratio) ENDDO END SUBROUTINE LMDZmoist_WRFmoist SUBROUTINE WRFmoist_LMDZmoist(wMOIST, dx, dy, dz, Nwmoist, wbdy, wqvid, wqcid, & & wqrid, Nlmdzq, lmdzMixRat) ! Subroutine to transform from WRF moisture array to the LMDZ moisture array IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy INTEGER, INTENT(IN) :: Nwmoist, wqvid, wqcid, wqrid INTEGER, INTENT(IN) :: Nlmdzq REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1,Nwmoist), INTENT(IN) :: wMOIST REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz,Nlmdzq), INTENT(OUT) :: lmdzMixRat ! Local INTEGER :: i, j, k, iq, dx2, dy2 INTEGER :: ddx, ddy ! INTEGER, ALLOCATABLE, DIMENSION(:) :: widlmdz ! REAL, DIMENSION(dx,dy) :: wcrMOIST REAL, DIMENSION(1-wbdy:dx-wbdy+1,1-wbdy:dy-wbdy+1) :: qcr_ratio, wcrMOIST !!!!!!! Variables ! wMOIST: WRF moisture array ! d[x,y,z]: Spatial dimension of the WRF moisture array ! Nwmoist: number of species in wMOIST ! wqvid: index for the water vapor mixing ratio ! wqcid: index for the cloud water mixing ratio ! wqrid: index for the rain water mixing ratio ! wcrMOIST: WRF array combination of cloud and rain liquid water ! Nlmdzq: number of species in LMDZ ! lmdxMixRat: LMDZ mixing ratio array ! widlmdz: vector with the id of the WRF tracers to be also included in the LMDZ moist array ddx=dx-2*wbdy ddy=dy-2*wbdy lmdzMixRat=0. DO k=1,dz ! Only taking water vapor and cloud water mixing ratios DO j=1,ddy DO i=1,ddx lmdzMixRat(ddx*(j-1)+i,k,1) = wMOIST(i,k,j,wqvid) lmdzMixRat(ddx*(j-1)+i,k,2) = wMOIST(i,k,j,wqcid) END DO END DO ENDDO !! IF (Nlmdzq > 2) THEN ! Whenever we would like to use more WRF water spicies some sort of mapping should be needed. ! Maybe something like this: ! IF (ALLOCATED(widlmdz)) DEALLOCATE(widlmdz) ! ALLOCATE(widlmdz(Nlmdzq-2)) !! DO iq=3,Nlmdzq !! DO k=1,dz !! lmdzMixRat(:,k,iq) = (iq-3.+(dz-k-1.)*1./(dz+1.))*1. !! END DO !! END DO !! END IF ! DEALLOCATE(widlmdz) END SUBROUTINE WRFmoist_LMDZmoist SUBROUTINE eta_to_pressure(etaP, psfc, ptop, dz, presP) ! Subroutine to transform from eta levels to pressure levels IMPLICIT NONE INTEGER, INTENT(IN) :: dz REAL, INTENT(IN) :: psfc, ptop REAL, DIMENSION(dz), INTENT(IN) :: etaP REAL, DIMENSION(dz), INTENT(OUT) :: presP !!!!!!Variables ! etaP: vector with the eta values ! psfc: pressure at the surface (in Pa) ! ptop: pressure at the top (in Pa) ! dz: number of vertical levels ! presP: equivalent vector of pressures (in Pa) presP = etaP*(psfc - ptop) + ptop END SUBROUTINE eta_to_pressure SUBROUTINE temp_to_theta(tempValues, pressValues, thetaValues, & & ids,ide, jds,jde, kds,kde, & & ims,ime, jms,jme, kms,kme, & & ips,ipe, jps,jpe, kps,kpe, & ! patch dims & i_start,i_end,j_start,j_end,kts,kte,num_tiles & ) ! Subroutine to transform from temperatures to WRF potential temperatures USE module_model_constants IMPLICIT NONE !-- thetaValues potential temperature values from WRF (300. to be added) !-- pressValues pressure !-- tempValues temperature values in K ! !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- ips start index for i in patch !-- ipe end index for i in patch !-- jps start index for j in patch !-- jpe end index for j in patch !-- kms start index for k in memory !-- kme end index for k in memory !-- i_start start indices for i in tile !-- i_end end indices for i in tile !-- j_start start indices for j in tile !-- j_end end indices for j in tile !-- kts start index for k in tile !-- kte end index for k in tile !-- num_tiles number of tiles ! !====================================================================== INTEGER, INTENT(IN) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & kts,kte, & num_tiles INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & & i_start,i_end,j_start,j_end REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: tempValues, & & pressValues REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: thetaValues ! Local REAL, PARAMETER :: pres0 = 100000. thetaValues = tempValues * ( pres0 / pressValues ) ** (r_d / cp) thetaValues = thetaValues - 300. END SUBROUTINE temp_to_theta SUBROUTINE theta_to_temp(thetaValues, pressValues, tempValues, & & ids,ide, jds,jde, kds,kde, & & ims,ime, jms,jme, kms,kme, & & ips,ipe, jps,jpe, kps,kpe, & ! patch dims & i_start,i_end,j_start,j_end,kts,kte,num_tiles & ) ! Subroutine to transform from WRF potential temperatures to temperature USE module_model_constants IMPLICIT NONE !-- thetaValues potential temperature values from WRF (300. to be added) !-- pressValues pressure !-- tempValues temperature values in K ! !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- ips start index for i in patch !-- ipe end index for i in patch !-- jps start index for j in patch !-- jpe end index for j in patch !-- kms start index for k in memory !-- kme end index for k in memory !-- i_start start indices for i in tile !-- i_end end indices for i in tile !-- j_start start indices for j in tile !-- j_end end indices for j in tile !-- kts start index for k in tile !-- kte end index for k in tile !-- num_tiles number of tiles ! !====================================================================== INTEGER, INTENT(IN) :: & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & kts,kte, & num_tiles INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & & i_start,i_end,j_start,j_end REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: thetaValues, & & pressvalues REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: tempValues ! Local REAL, PARAMETER :: pres0 = 100000. INTEGER :: i,j,k !! PRINT *,' thetaValues: ', LBOUND(thetaValues),' x ', UBOUND(thetaValues) !! DO k=kms,kme-2 !! DO j=jms,jme !! DO i=ims,ime !! IF (thetaValues(i,k,j) == 0. .AND. i > 0 .AND. i < 32 & !! .AND. j > 0 .AND. j < 32) THEN !! PRINT *,' Zth: T ',i,k,j,thetaValues(i,k,j) !! END IF !! END DO !! END DO !! END DO tempValues = ( thetaValues + 300.)*( pressValues / pres0 ) ** (r_d / cp) !! DO k=kms,kme-2 END SUBROUTINE theta_to_temp SUBROUTINE temp_to_theta1(tempValue, pressValue, thetaValue) ! Subroutine to transform from temperature to potential temperature USE module_model_constants IMPLICIT NONE !!!!!!! Variables ! thetaValue potential temperature value ! pressValue pressure ! tempValue temperature value in K ! REAL, INTENT(IN) :: tempValue, pressValue REAL, INTENT(OUT) :: thetaValue ! Local REAL, PARAMETER :: pres0 = 100000. thetaValue = tempValue * ( pres0 / pressValue ) ** (r_d / cp) END SUBROUTINE temp_to_theta1 SUBROUTINE theta_to_temp1(thetaValue, pressValue, tempValue) ! Subroutine to transform from potential temperature to temperature USE module_model_constants IMPLICIT NONE !!!!!!! Variables ! thetaValue potential temperature ! pressValue pressure ! tempValue temperature value in K REAL, INTENT(IN) :: thetaValue, pressvalue REAL, INTENT(OUT) :: tempValue ! Local REAL, PARAMETER :: pres0 = 100000. tempValue = thetaValue*( pressValue / pres0 ) ** (r_d / cp) END SUBROUTINE theta_to_temp1 SUBROUTINE def_get_scalar_value(filen, valn, ds, di, dr, db) ! Subroutine to get a scalar value from a .def file ! NOTE: For that scalar values with scientific notation (EN, or e-N) it will not work... ! next version! IMPLICIT NONE CHARACTER(LEN=100), INTENT(IN) :: filen, valn CHARACTER(LEN=100), INTENT(OUT) :: ds INTEGER, INTENT(OUT) :: di REAL, INTENT(OUT) :: dr LOGICAL, INTENT(OUT) :: db ! Local INTEGER :: funit, ios, Lvaln INTEGER :: iline, lineval CHARACTER(LEN=1) :: kindvalue CHARACTER(LEN=100) :: varname,varname0 CHARACTER(LEN=50) :: subname, errmsg, value, & value0 LOGICAL :: fexists, is_used !!!!!!! Variables ! filen: Name of the def file to get the value ! valn: name of the value to get ! ds: string value (up to 100 characters) ! di: integer value ! dr: real value ! db: logical value ! value: value get from the file ! kindvalue: kind of the value ('s': string, 'i': integer, 'r': real, 'b':boolean ) subname = 'def_get_scalar_value' errmsg = 'ERROR -- error -- ERROR -- error' INQUIRE(FILE=TRIM(filen), EXIST=fexists) IF (fexists) THEN Lvaln = LEN_TRIM(valn) ds = 'None' di = 0 dr = 0. db = .FALSE. DO funit=10,100 INQUIRE(UNIT=funit, OPENED=is_used) IF (.NOT. is_used) EXIT END DO OPEN(funit, FILE=TRIM(filen), STATUS='old', FORM='formatted', IOSTAT=ios) IF ( ios /= 0 ) THEN PRINT *,errmsg PRINT *, ' ' // TRIM(subname) // ": ERROR opening '" // TRIM(filen) // "'" PRINT *,'ios: ',ios STOP END IF ! iunit = get_unused_unit() ! IF ( iunit <= 0 ) THEN ! IF ( wrf_dm_on_monitor() ) THEN ! CALL wrf_error_fatal('Error in module_ra_cam_support: could not find a free Fortran unit.') ! END IF ! END IF lineval = 0 DO iline=1, 10000 READ(funit,*,END=125)varname0 varname = TRIM(varname0) IF (varname(1:Lvaln) == TRIM(valn)) THEN lineval = iline EXIT END IF END DO 125 CONTINUE IF (lineval == 0) THEN PRINT *,errmsg PRINT *,' ' // TRIM(subname) // ": In file '" // TRIM(filen) // & & "' does not exist a variable called '" // TRIM(valn) // '" !!!' STOP ELSE CLOSE(funit) OPEN(funit, FILE=TRIM(filen), STATUS='old', FORM='formatted', IOSTAT=ios) DO iline=1,lineval-1 READ(funit,*)varname END DO READ(funit,*)value0 CLOSE(UNIT=funit) CALL string_split(value0, '=', 2, value) CALL string_split(value, ':', 2, value) CALL value_kind(value, kindvalue) ! PRINT *,'value: ',TRIM(value), ' kind: ',kindvalue IF (kindvalue == 'i') THEN READ(value, '(I50)')di ELSE IF (kindvalue == 'r') THEN READ(value, '(F50.25)')dr ELSE IF (kindvalue == 'e') THEN READ (value, '(E50.25)')dr ELSE IF (kindvalue == 's') THEN ds = TRIM(value) ELSE IF (TRIM(value) == 'y' .OR. TRIM(value) == 'yes' .OR. TRIM(value) == 'Y' & .OR. TRIM(value) == 'Yes' .OR. TRIM(value) == 'YES') db = .TRUE. IF (TRIM(value) == 'n' .OR. TRIM(value) == 'no' .OR. TRIM(value) == 'N' & .OR. TRIM(value) == 'No' .OR. TRIM(value) == 'NO') db = .FALSE. END IF END IF ELSE PRINT *,errmsg PRINT *,' ' // TRIM(subname) // ": File '" // TRIM(filen) // "' does not exist !!" STOP END IF END SUBROUTINE def_get_scalar_value SUBROUTINE value_kind(value, kindval) ! Subroutine to determine which kind of value is given IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: value CHARACTER(LEN=1), INTENT(OUT) :: kindval ! Local INTEGER :: ic INTEGER :: Lval, Nc, NnumVal, NstrVal LOGICAL :: is_point, is_exp !!!!!!! Variables ! value: value to determine the cahracteristics ! kindval: kind of value. 's': string, 'i': integer, 'r': real, 'e': real scientific notation, ! 'b':boolean ! Lval: length of characters of the value ! NnumVal: number of numeric characers ! NstrVal: number of string characters ! is_point: boolean variable indicating the presence of a '.' ! is_exp: boolean variable indicating the presence of a 'E' or 'e' (exponential) IF (TRIM(value) == 'y' .OR. TRIM(value) == 'n' .OR. TRIM(value) == 'yes' .OR. & TRIM(value) == 'no' .OR. TRIM(value) == 'Y' .OR. TRIM(value) == 'N' .OR. & TRIM(value) == 'Yes' .OR. TRIM(value) == 'No' .OR. TRIM(value) == 'YES' .OR. & TRIM(value) == 'NO') THEN ! PRINT *,' boolean!' kindval = 'b' ELSE Lval = LEN_TRIM(value) NnumVal = 0 NstrVal = 0 is_point = .FALSE. DO ic=1,Lval Nc = ICHAR(value(ic:ic)) ! PRINT *, 'char: ',value(ic:ic), ' Nc:',Nc IF (Nc >= 48 .AND. Nc <= 57) THEN Nnumval = Nnumval + 1 ELSE NstrVal = NstrVal + 1 IF (value(ic:ic) == '.') is_point = .TRUE. IF ((value(ic:ic) == 'e') .OR. (value(ic:ic) == 'E')) is_exp = .TRUE. END IF END DO ! PRINT *,'NnumVal:', NnumVal, 'NstrVal: ', NstrVal, 'point: ', is_point IF (Nnumval == Lval) THEN ! PRINT *,' integer!' kindval = 'i' ELSE IF ( (Nnumval == Lval -1) .AND. (is_point .EQV. .TRUE.) .AND. & (NstrVal == 1)) THEN ! PRINT *,' real!' kindval = 'r' ELSE IF ( ((Nnumval == Lval -2) .OR.(Nnumval == Lval -3)) .AND. & (is_point .EQV. .TRUE.) .AND. (is_exp .EQV. .TRUE.) .AND. & ((NstrVal == 2) .OR. (NstrVal == 3)) ) THEN ! PRINT *,' real exponential!' kindval = 'e' ELSE ! PRINT *,' string!' kindval = 's' END IF END IF END SUBROUTINE value_kind SUBROUTINE string_split(string, charsplit, Nsec, sec) ! Subtoutine to split a string according to a character and returns the 'Nsec' of the split IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: string CHARACTER(LEN=1), INTENT(IN) :: charsplit INTEGER, INTENT(IN) :: Nsec CHARACTER(LEN=50), INTENT(OUT) :: sec ! Local INTEGER :: ic INTEGER :: poschar, Lstring INTEGER :: isec, inisec CHARACTER(LEN=50) :: str, secchar !!!!!!! Variables ! string: string to split ! charsplit: character to use to split ! Nsec: number of section to return ! sec: returned section of the string ! PRINT *,'string: ',string str = TRIM(string) poschar=0 Lstring = LEN_TRIM(str) isec = 1 poschar = INDEX(str, charsplit) inisec = 1 IF (poschar == 0) THEN sec = str ELSE DO ic=poschar,50 secchar = str(inisec:poschar-1) IF (isec == Nsec) sec = secchar inisec = poschar + 1 poschar = INDEX(str(inisec:Lstring), charsplit) IF (poschar == 0) poschar = Lstring + 2 isec = isec + 1 END DO END IF END SUBROUTINE string_split SUBROUTINE get_lmdz(ddx,ddy,ddz,bdy,dlmdz,Nks,NzO,plon,plat,pzmasq,ppctsrf,pftsol, & pftsoil,pqsol,pqsurf,pfalb1,pfalb2,pevap,psnow,pradsol,psolsw,psollw,pfder, & prain_fall,psnow_fall,pagesno,pzmea,pzstd,pzgam,pzthe,pzpic,pzval,prugoro, & prugos,pzmax0,pf0,pwake_s,pwake_cstar,pwake_pe,pwake_fip,psst, & palbedo,pt_ancien,pq_ancien,pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs, & pema_work1,pema_work2,pwake_deltat,pwake_deltaq,pfm_therm,pentr_therm, & pdetr_therm,pnat,pwrf_grid) ! Subroutine to get LMDZ variables with values from WRF ! WRF modules USE module_domain_type ! LMDZ modules USE indice_sol_mod IMPLICIT NONE TYPE(domain) , TARGET :: pwrf_grid INTEGER, INTENT(IN) :: ddx, ddy, ddz, bdy, & dlmdz, Nks, NzO REAL, DIMENSION(dlmdz), INTENT(IN) :: plon, plat, pzmasq, & pqsol, pradsol, psolsw, psollw, pfder, prain_fall, psnow_fall, pzmea, pzstd, & pzgam, pzthe, pzpic, pzval, prugoro, pzmax0, pf0, pwake_s, & pwake_cstar, pwake_pe, pwake_fip, psst, palbedo, pnat ! What to do with these variables? ,pnat,pbils,prug REAL, DIMENSION(dlmdz,Nks), INTENT(IN) :: ppctsrf,pftsol,pqsurf, & pfalb1,pfalb2,pevap,psnow,prugos,pagesno REAL, DIMENSION(dlmdz,ddz), INTENT(IN) :: pt_ancien,pq_ancien, & pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,pema_work1,pema_work2,pwake_deltat,& pwake_deltaq,pentr_therm,pdetr_therm REAL, DIMENSION(dlmdz,NzO,Nks), INTENT(IN) :: pftsoil REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN) :: pfm_therm ! Local INTEGER :: ix,iy,iz,ixy,dx,dy CHARACTER(LEN=50) :: LMDZvarmethod ! Filling all 2D variables !! DO iy=1,ddy DO ix=1,ddx ixy=ddx*(iy-1)+ix ! It does not change on time: pzmask ! It might change due to ice dynamics? pwrf_grid%lter(ix,iy) = ppctsrf(ixy,is_ter) pwrf_grid%llic(ix,iy) = ppctsrf(ixy,is_lic) pwrf_grid%loce(ix,iy) = ppctsrf(ixy,is_oce) pwrf_grid%lsic(ix,iy) = ppctsrf(ixy,is_sic) DO iz=1,Nks pwrf_grid%ltksoil(ix,iz,iy) = pftsol(ixy,iz) END DO DO iz=1,NzO pwrf_grid%lotter(ix,iz,iy) = pftsoil(ixy,iz,1) pwrf_grid%lotlic(ix,iz,iy) = pftsoil(ixy,iz,2) pwrf_grid%lotoce(ix,iz,iy) = pftsoil(ixy,iz,3) pwrf_grid%lotsic(ix,iz,iy) = pftsoil(ixy,iz,4) END DO DO iz=1,Nks pwrf_grid%lqksoil(ix,iz,iy) = pqsurf(ixy,iz) END DO pwrf_grid%lwsol(ix,iy) = pqsol(ixy) DO iz=1,Nks pwrf_grid%lalbksoil(ix,iz,iy) = pfalb1(ixy,iz) pwrf_grid%llwalbksoil(ix,iz,iy) = pfalb2(ixy,iz) pwrf_grid%levapksoil(ix,iz,iy) = pevap(ixy,iz) pwrf_grid%lsnowksoil(ix,iz,iy) = psnow(ixy,iz) END DO pwrf_grid%lrads(ix,iy) = pradsol(ixy) pwrf_grid%lsolsw(ix,iy) = psolsw(ixy) pwrf_grid%lsollw(ix,iy) = psollw(ixy) pwrf_grid%lfder(ix,iy) = pfder(ixy) pwrf_grid%lrain(ix,iy) = prain_fall(ixy) pwrf_grid%lsnow(ix,iy) = psnow_fall(ixy) DO iz=1,Nks pwrf_grid%lrugksoil(ix,iz,iy) = prugos(ixy,iz) pwrf_grid%lagesnoksoil(ix,iz,iy) = pagesno(ixy,iz) END DO ! These variables do not change on time! !! ! pzmea, pzstd, pzsig, pzgam, pzthe, pzpic, pzval, prugsrel pwrf_grid%lrugsea(ix,iy) = prugos(ixy,is_oce) ! pwrf_grid%lrunofflic(ix,iy)= prun_off_lic(ixy) pwrf_grid%lzmax0(ix,iy) = pzmax0(ixy) pwrf_grid%lf0(ix,iy) = pf0(ixy) pwrf_grid%lwake_s(ix,iy) = pwake_s(ixy) pwrf_grid%lwake_cstar(ix,iy) = pwake_cstar(ixy) pwrf_grid%lwake_pe(ix,iy) = pwake_pe(ixy) pwrf_grid%lwake_fip(ix,iy) = pwake_fip(ixy) pwrf_grid%lnat(ix,iy) = pnat(ixy) pwrf_grid%sst(ix,iy) = psst(ixy) ! pwrf_grid%lbils(ix,iy) = pbils(ixy) pwrf_grid%albedo(ix,iy) = palbedo(ixy) ! pwrf_grid%lrug(ix,iy) = prug(ixy) END DO END DO dx=ddx+2*bdy dy=ddy+2*bdy CALL vect_mat(pclwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon) CALL vect_mat(prnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon) CALL vect_mat(pratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs) CALL vect_mat(pema_work1, dx, dy, ddz, bdy, pwrf_grid%lema_work1) CALL vect_mat(pema_work2, dx, dy, ddz, bdy, pwrf_grid%lema_work2) CALL vect_mat(pwake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat) CALL vect_mat(pwake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq) CALL vect_mat(pfm_therm, dx, dy, ddz+1, bdy, pwrf_grid%lfm_therm) CALL vect_mat(pentr_therm, dx, dy, ddz, bdy, pwrf_grid%lentr_therm) CALL vect_mat(pdetr_therm, dx, dy, ddz, bdy, pwrf_grid%ldetr_therm) RETURN END SUBROUTINE get_lmdz SUBROUTINE get_lowbdy(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, pkglo, pwrf_grid, & ppctsrf, pftsol, prugos, palbedo, psst, ptsurf_lim, pz0_lim, palb_lim, prug_glo, & ppct_glo) ! Subroutine to load LMDZ limit variables with values from lowbdy WRF ! WRF modules ! USE module_model_constants USE module_domain_type ! LMDZ modules USE indice_sol_mod IMPLICIT NONE TYPE(domain) , TARGET :: pwrf_grid INTEGER, INTENT(IN) :: ddx, ddy, ddz, bdy, & dlmdz, Nks, Nz0, pkglo REAL, DIMENSION(dlmdz), INTENT(INOUT) :: psst, palbedo, & ptsurf_lim, pz0_lim, palb_lim ! What to do with these variables? ,pnat,pbils REAL, DIMENSION(dlmdz,Nks), INTENT(INOUT) :: ppctsrf, pftsol, prugos REAL, DIMENSION(pkglo), INTENT(OUT) :: prug_glo REAL, DIMENSION(pkglo,Nks), INTENT(OUT) :: ppct_glo ! Local INTEGER :: ix,iy,iz,ixy,dx,dy,ixx,iyy CHARACTER(LEN=50) :: LMDZvarmethod,fname,errmsg fname='get_lowbdy' errmsg='ERROR -- error -- ERROR -- error' ! Filling all 2D variables !! DO iy=1,ddy DO ix=1,ddx ixy=ddx*(iy-1)+ix ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy) ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy) IF ( pwrf_grid%xice(ix,iy) /= pwrf_grid%lsic(ix,iy) ) THEN ! No ocean/frozen ocean fraction before IF (pwrf_grid%loce(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,3,iy) = 0.0001 IF (pwrf_grid%lsic(ix,iy) == 0.) pwrf_grid%lrugksoil(ix,4,iy) = 0.001 pwrf_grid%loce(ix,iy) = 1. - pwrf_grid%lsic(ix,iy) pwrf_grid%lsic(ix,iy) = pwrf_grid%lsic(ix,iy) ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy) ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy) ! Recomputing total roughness !! pwrf_grid%lrug(ix,iy)=pwrf_grid%lter(ix,iy)*pwrf_grid%lrugksoil(ix,1,iy) +& pwrf_grid%llic(ix,iy)*pwrf_grid%lrugksoil(ix,2,iy) + & pwrf_grid%loce(ix,iy)*pwrf_grid%lrugksoil(ix,3,iy) + & pwrf_grid%lsic(ix,iy)*pwrf_grid%lrugksoil(ix,4,iy) END IF ! Full SST grid points must have a SST value from the wrfinput (WRF 1/0 landmask) ! but that LMDZ points with fractions of land/mask could not. Using the closest one ! if it is an isolated grid point, using TSK instead IF (pwrf_grid%sst(ix,iy) < 200.15) THEN pwrf_grid%ltksoil(ix,3,iy) = -9. IF ( (ix > 1).AND.(ix < ddx) .AND. (iy > 1).AND.(iy < ddy) ) THEN DO ixx=-1,1 DO iyy=-1,1 ! PRINT *,ixx,iyy,wrf_grid%sst(ix+ixx,iy+iyy) IF (pwrf_grid%sst(ix+ixx,iy+iyy) > 200.15) THEN pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix+ixx,iy+iyy) EXIT END IF END DO IF ( pwrf_grid%ltksoil(ix,3,iy) > 200.15) EXIT END DO IF ( pwrf_grid%ltksoil(ix,3,iy) == -9.) pwrf_grid%ltksoil(ix,3,iy) = & pwrf_grid%tsk(ix,iy) ELSE pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%tsk(ix,iy) END IF ELSE pwrf_grid%ltksoil(ix,3,iy) = pwrf_grid%sst(ix,iy) END IF IF ( (pwrf_grid%ltksoil(ix,3,iy) < 200.15) .OR. & (pwrf_grid%ltksoil(ix,3,iy) > 370.) ) THEN PRINT *,TRIM(errmsg) ! WRITE(wrf_err_message,*) ' ' // TRIM(fname) // & ! ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ', & ! wrf_grid%loce(ix,iy),' has a tsoil(oce)= ', & ! wrf_grid%ltksoil(ix,3,iy),' sst: ',wrf_grid%sst(ix,iy),' tsk: ', & ! wrf_grid%tsk(ix,iy) ! CALL wrf_error_fatal(TRIM(wrf_err_message)) PRINT *, ' ' // TRIM(fname) // & ' SST inconsistency at point: ', ix, ', ', iy,' with loce: ', & pwrf_grid%loce(ix,iy),' has a tsoil(oce)= ', & pwrf_grid%ltksoil(ix,3,iy),' sst: ', pwrf_grid%sst(ix,iy),' tsk: ', & pwrf_grid%tsk(ix,iy) END IF psst(ixy) = pwrf_grid%ltksoil(ix,3,iy) palbedo(ixy) = pwrf_grid%albbck(ix,iy) DO iz=1,Nks prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy) END DO ! In ocean_forced_mod is used as tsurf_lim = sst(knindex(1:knon)). Let's complicate it ! a little bit... ptsurf_lim(ixy)= psst(ixy) ! In surf_land_bucket_mod is used as z0_new=rugos(knindex(1:knon),1) [z0_new == z0_limit] pz0_lim(ixy) = prugos(ixy,1) palb_lim(ixy) = palbedo(ixy) prug_glo(ixy) = prugos(ixy,is_ter) ppct_glo(ixy,:) = ppctsrf(ixy,:) END DO END DO RETURN END SUBROUTINE get_lowbdy SUBROUTINE load_lmdz(pwrf_grid,ddx,ddy,ddz,bdy,dlmdz,Nks,NzO,prlon,prlat,pzmasq, & ppctsrf,pftsol,pftsoil,pqsurf,pqsol,pfalb1,pfalb2,pevap,psnow,pradsol,psolsw, & psollw,pfder,prain_fall,psnow_fall,pagesno,pzmea,pzstd,pzgam,pzthe,pzpic,pzval, & prugoro,prugos,prun_off_lic,pzmax0,pf0,pwake_s,pwake_cstar,pwake_pe,pwake_fip, & psst,palbedo,pt_ancien,pq_ancien,pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs, & pema_work1,pema_work2,pwake_deltat,pwake_deltaq,pfm_therm,pentr_therm, & pdetr_therm,pnat) ! Subroutine to load LMDZ variables with values from WRF ! WRF modules ! USE module_model_constants USE module_domain_type ! LMDZ modules USE indice_sol_mod IMPLICIT NONE TYPE(domain) , TARGET :: pwrf_grid INTEGER, INTENT(IN) :: ddx, ddy, ddz, bdy, & dlmdz, Nks, NzO REAL, DIMENSION(dlmdz), INTENT(OUT) :: prlon, prlat, pzmasq, & pqsol, pradsol, psolsw, psollw, pfder, prain_fall, psnow_fall, pzmea, pzstd, & pzgam, pzthe, pzpic, pzval, prugoro, prun_off_lic, pzmax0, pf0, pwake_s, & pwake_cstar, pwake_pe, pwake_fip, psst, palbedo,pnat REAL, DIMENSION(dlmdz,Nks), INTENT(OUT) :: ppctsrf,pftsol,pqsurf, & pfalb1,pfalb2,pevap,psnow,prugos,pagesno REAL, DIMENSION(dlmdz,ddz), INTENT(OUT) :: pt_ancien,pq_ancien, & pu_ancien,pv_ancien,pclwcon,prnebcon,pratqs,pema_work1,pema_work2,pwake_deltat,& pwake_deltaq,pentr_therm,pdetr_therm REAL, DIMENSION(dlmdz,NzO,Nks), INTENT(OUT) :: pftsoil REAL, DIMENSION(dlmdz,ddz+1), INTENT(OUT) :: pfm_therm ! Local INTEGER :: ix,iy,iz,ixy,dx,dy CHARACTER(LEN=50) :: LMDZvarmethod,fname fname='load_lmdz' ! Filling all 2D variables !! PRINT *,' ' // TRIM(fname) // '; LUBOUND(grid%lmsk): ',LBOUND(pwrf_grid%lmsk), UBOUND(pwrf_grid%lmsk) PRINT *,' ' // TRIM(fname) // '; ddx ddy: ',ddy,ddx DO iy=1,ddy DO ix=1,ddx ixy=ddx*(iy-1)+ix prlon(ixy) = pwrf_grid%xlong(ix,iy) prlat(ixy) = pwrf_grid%xlat(ix,iy) pzmasq(ixy) = pwrf_grid%lmsk(ix,iy) ppctsrf(ixy,is_ter) = pwrf_grid%lter(ix,iy) ppctsrf(ixy,is_lic) = pwrf_grid%llic(ix,iy) ppctsrf(ixy,is_oce) = pwrf_grid%loce(ix,iy) ppctsrf(ixy,is_sic) = pwrf_grid%lsic(ix,iy) DO iz=1,Nks pftsol(ixy,iz) = pwrf_grid%ltksoil(ix,iz,iy) END DO DO iz=1,NzO pftsoil(ixy,iz,1) = pwrf_grid%lotter(ix,iz,iy) pftsoil(ixy,iz,2) = pwrf_grid%lotlic(ix,iz,iy) pftsoil(ixy,iz,3) = pwrf_grid%lotoce(ix,iz,iy) pftsoil(ixy,iz,4) = pwrf_grid%lotsic(ix,iz,iy) END DO DO iz=1,Nks pqsurf(ixy,iz) = pwrf_grid%lqksoil(ix,iz,iy) END DO pqsol(ixy) = pwrf_grid%lwsol(ix,iy) DO iz=1,Nks pfalb1(ixy,iz) = pwrf_grid%lalbksoil(ix,iz,iy) pfalb2(ixy,iz) = pwrf_grid%llwalbksoil(ix,iz,iy) pevap(ixy,iz) = pwrf_grid%levapksoil(ix,iz,iy) psnow(ixy,iz) = pwrf_grid%lsnowksoil(ix,iz,iy) END DO ! No way to get back the independent radiative values (unless using LMDZ 'hist*' ! files). Use it in the LMDZ's added way also from grid (no need to initialize it, ! no radiation on t=0) pradsol(ixy) = pwrf_grid%lrads(ix,iy) psolsw(ixy) = pwrf_grid%lsolsw(ix,iy) psollw(ixy) = pwrf_grid%lsollw(ix,iy) pfder(ixy) = pwrf_grid%lfder(ix,iy) ! LMDZ does not provide a way (in the restart) to distingusih between convective, ! non-convective. Thus Use it as a whole prain_fall(ixy) = pwrf_grid%lrain(ix,iy) psnow_fall(ixy) = pwrf_grid%lsnow(ix,iy) DO iz=1,Nks prugos(ixy,iz) = pwrf_grid%lrugksoil(ix,iz,iy) pagesno(ixy,iz) = pwrf_grid%lagesnoksoil(ix,iz,iy) END DO pzmea(ixy) = pwrf_grid%lzmea(ix,iy) pzstd(ixy) = pwrf_grid%lzstd(ix,iy) pzgam(ixy) = pwrf_grid%lzgam(ix,iy) pzthe(ixy) = pwrf_grid%lzthe(ix,iy) pzpic(ixy) = pwrf_grid%lzpic(ix,iy) pzval(ixy) = pwrf_grid%lzval(ix,iy) prugoro(ixy) = pwrf_grid%lzrugsrel(ix,iy) ! prugos(ixy,is_oce) = pwrf_grid%lrugsea(ix,iy) prun_off_lic(ixy) = pwrf_grid%lrunofflic(ix,iy) pzmax0(ixy) = pwrf_grid%lzmax0(ix,iy) pf0(ixy) = pwrf_grid%lf0(ix,iy) pwake_s(ixy) = pwrf_grid%lwake_s(ix,iy) pwake_cstar(ixy) = pwrf_grid%lwake_cstar(ix,iy) pwake_pe(ixy) = pwrf_grid%lwake_pe(ix,iy) pwake_fip(ixy) = pwrf_grid%lwake_fip(ix,iy) !! limit.nc ! nat == pctsrf, not needed ? pnat(ixy) = pwrf_grid%lnat(ix,iy) psst(ixy) = pwrf_grid%sst(ix,iy) ! not used, not defined.... ! pbils(ixy) = pwrf_grid%lbils(ix,iy) palbedo(ixy) = pwrf_grid%albedo(ix,iy) ! already done in previous steps ! prugos(ixy) = pwrf_grid%lrug(ix,iy) END DO END DO dx=ddx+2*bdy dy=ddy+2*bdy CALL mat_vect(pwrf_grid%t_2, dx, dy, ddz, bdy, pt_ancien) ! New variable added in the Registry qv_2 CALL mat_vect(pwrf_grid%qv_2, dx, dy, ddz, bdy, pq_ancien) CALL mat_vect(pwrf_grid%u_2, dx, dy, ddz, bdy, pu_ancien) CALL mat_vect(pwrf_grid%v_2, dx, dy, ddz, bdy, pv_ancien) CALL mat_vect(pwrf_grid%lclwcon, dx, dy, ddz, bdy, pclwcon) CALL mat_vect(pwrf_grid%lrnebcon, dx, dy, ddz, bdy, prnebcon) CALL mat_vect(pwrf_grid%lratqs, dx, dy, ddz, bdy, pratqs) CALL mat_vect(pwrf_grid%lema_work1, dx, dy, ddz, bdy, pema_work1) CALL mat_vect(pwrf_grid%lema_work2, dx, dy, ddz, bdy, pema_work2) CALL mat_vect(pwrf_grid%lwake_deltat, dx, dy, ddz, bdy, pwake_deltat) CALL mat_vect(pwrf_grid%lwake_deltaq, dx, dy, ddz, bdy, pwake_deltaq) CALL mat_vect(pwrf_grid%lfm_therm, dx, dy, ddz+1, bdy, pfm_therm) CALL mat_vect(pwrf_grid%lentr_therm, dx, dy, ddz, bdy, pentr_therm) CALL mat_vect(pwrf_grid%ldetr_therm, dx, dy, ddz, bdy, pdetr_therm) RETURN END SUBROUTINE load_lmdz SUBROUTINE vect_mat(vect, dx, dy, dz, wbdy, mat) ! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat ! Local INTEGER :: i,j,k INTEGER :: ddx, ddy INTEGER :: lp, il, jl ! Variables ! d[x/y/z]: dimensions ! mat: matrix to transform ! vect: vector to produce mat=0. ddx=dx-2*wbdy ddy=dy-2*wbdy ! lp=321 ! jl=INT(lp/ddx)+1 ! il=lp-ddx*(jl-1) DO k=1,dz DO j=1,ddy DO i=1,ddx mat(i,k,j)=vect(ddx*(j-1)+i,k) END DO END DO END DO END SUBROUTINE vect_mat SUBROUTINE vect_mat_zstagg(vect, dx, dy, dz, wbdy, mat) ! Subroutine to transform a LMDZ physic 1D vector to a 3D matrix IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(IN) :: vect REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(OUT) :: mat ! Local INTEGER :: i,j,k INTEGER :: ddx, ddy ! Variables ! d[x/y/z]: dimensions ! mat: matrix to transform ! vect: vector to produce ddx=dx-2*wbdy ddy=dy-2*wbdy mat=0. DO k=1,dz DO j=1,ddy DO i=1,ddx mat(i,k,j)=vect(ddx*(j-1)+i,k) END DO END DO END DO END SUBROUTINE vect_mat_zstagg SUBROUTINE mat_vect(mat, dx, dy, dz, wbdy, vect) ! Subroutine to transform a 3D matrix to a LMDZ physic 1D vector IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz+1,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) :: & vect ! Local INTEGER :: i,j,k INTEGER :: ddx,ddy ! Variables ! d[x/y/z]: dimensions ! mat: matrix to transform ! vect: vector to produce ! Removing HALO ddx=dx-2*wbdy ddy=dy-2*wbdy vect=0. DO k=1,dz DO j=1,ddy DO i=1,ddx vect(ddx*(j-1)+i,k)= mat(i,k,j) END DO END DO END DO END SUBROUTINE mat_vect SUBROUTINE mat_vect_zstagg(mat, dx, dy, dz, wbdy, vect) ! Subroutine to transform a 3D matrix z-staggered to a LMDZ physic 1D vector IMPLICIT NONE INTEGER, INTENT(IN) :: dx, dy, dz, wbdy REAL, DIMENSION(1-wbdy:dx-wbdy+1,dz,1-wbdy:dy-wbdy+1), INTENT(IN) :: mat REAL, DIMENSION((dx-2*wbdy)*(dy-2*wbdy),dz), INTENT(OUT) :: & vect ! Local INTEGER :: i,j,k INTEGER :: ddx,ddy ! Variables ! d[x/y/z]: dimensions ! mat: matrix to transform ! vect: vector to produce ! Removing HALO ddx=dx-2*wbdy ddy=dy-2*wbdy DO k=1,dz DO j=1,ddy DO i=1,ddx vect(ddx*(j-1)+i,k)= mat(i,k,j) END DO END DO END DO END SUBROUTINE mat_vect_zstagg SUBROUTINE WRFlanduse_LMDZKsoil(lndcn,is,ie,js,je,dx,dy,ncat,lndu,mask,kt,kl,ko,ks) ! Subroutine to retrieve LMDZ Ksoil classification using WRF landuse IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: lndcn INTEGER, INTENT(IN) :: is,ie,js,je,dx,dy,ncat REAL, DIMENSION(is:ie,ncat,js:je), INTENT(IN) :: lndu REAL, DIMENSION(is:ie,js:je), INTENT(IN) :: mask REAL, DIMENSION(is:ie,js:je), INTENT(OUT) :: kt,kl,ko,ks ! Local INTEGER :: i,j,k,ij,Nnones,Ncattot CHARACTER(LEN=50) :: fname, errmsg, wrnmsg REAL, DIMENSION(is:ie,js:je) :: totlndu INTEGER, ALLOCATABLE, DIMENSION(:) :: ter_wk, lic_wk, oce_wk, & sic_wk LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: wk INTEGER :: Nter_wk, Nlic_wk, & Noce_wk, Nsic_wk, klks_equal REAL*8, DIMENSION(is:ie,js:je) :: Ksoilval !!!!!!!! Variables ! lndcn: land-use data-set (same as in VEGPARM.TBL) ! dx,dy: horizontal dimension of the domain ! ncat: number of categories of the data-set ! lndu: WRF landuse matrix ! mask: land-sea mask from WRF ! k[t/l/o/s]: LMDZ ter, lic, oce, sic matrices ! Ncattot: theoretical number of laduses for a given data-base ! [ter/lic/oce/sic]_wk: vectors with the indices of the landuse for LMDZ Ksoil ! N[ter/lic/oce/sic]_wk: number of indices of the landuse for LMDZ Ksoil ! wk: boolean matrix for each Ksoil indicating 'true' if the WRF landuse is for Ksoil fname="'WRFlanduse_LMDZKsoil'" errmsg='ERROR -- error -- ERROR -- error' wrnmsg='WARNING -- warning -- WARNING -- warning' ! Checking if at all the grid points SUM(landuse) == 1. !! PRINT *,' ' // TRIM(fname) // " using '" // TRIM(lndcn) // "' WRF land types" totlndu = SUM(lndu, 2) Nnones=0 DO i=1,dx DO j=1,dy IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1 END DO END DO IF (Nnones /= 0) THEN PRINT *,TRIM(errmsg) PRINT *,' ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' // & 'TOTlanduse == 1. from WRF !!!!!' DO i=1,dx DO j=1,dy IF (totlndu(i,j) /= 1.) PRINT *,' ',i,', ',j,' total land use: ', & totlndu(i,j),' indiv. values :', lndu(i,:,j) END DO END DO END IF SELECT CASE (lndcn) CASE ('OLD') Ncattot = 13 Nter_wk = 11 Nlic_wk = 1 Noce_wk = 1 Nsic_wk = 1 IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk) ALLOCATE(ter_wk(Nter_wk)) IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk) ALLOCATE(lic_wk(Nlic_wk)) IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk) ALLOCATE(oce_wk(Noce_wk)) IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk) ALLOCATE(sic_wk(Nsic_wk)) ! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3) !! ! 1,'Urban land' 2,'Agriculture' ! 3,'Range-grassland' 4,'Deciduous forest' ! 5,'Coniferous forest' 6,'Mixed forest and wet land' ! 7,'Water' 8,'Marsh or wet land' ! 9,'Desert' 10,'Tundra' ! 11,'Permanent ice' 12,'Tropical or subtropical forest' ! 13,'Savannah' ter_wk = (/ 1,2,3,4,5,6,8,9,10,12,13 /) lic_wk = (/ 11 /) oce_wk = (/ 7 /) sic_wk = (/ 11 /) CASE ('USGS') Ncattot = 33 Nter_wk = 22 Nlic_wk = 1 Noce_wk = 1 Nsic_wk = 1 IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk) ALLOCATE(ter_wk(Nter_wk)) IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk) ALLOCATE(lic_wk(Nlic_wk)) IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk) ALLOCATE(oce_wk(Noce_wk)) IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk) ALLOCATE(sic_wk(Nsic_wk)) ! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3) !! ! 1,'Urban and Built-Up Land' 2,'Dryland Cropland and Pasture' ! 3,'Irrigated Cropland and Pasture' 4,'Mixed Dryland/Irrigated Cropland and Pasture' ! 5,'Cropland/Grassland Mosaic' 6,'Cropland/Woodland Mosaic' ! 7,'Grassland' 8,'Shrubland' ! 9,'Mixed Shrubland/Grassland' 10,'Savanna' !11,'Deciduous Broadleaf Forest' 12,'Deciduous Needleleaf Forest' !13,'Evergreen Broadleaf Forest' 14,'Evergreen Needleleaf Forest' !15,'Mixed Forest' 16,'Water Bodies' !17,'Herbaceous Wetland' 18,'Wooded Wetland' !19,'Barren or Sparsely Vegetated' 20,'Herbaceous Tundra' !21,'Wooded Tundra' 22,'Mixed Tundra' !23,'Bare Ground Tundra' 24,'Snow or Ice' !25,'Playa' 26,'Lava' !27,'White Sand' 28,'Unassigned' !29,'Unassigned' 30,'Unassigned' !31,'Low Intensity Residential ' 32,'High Intensity Residential' !33,'Industrial or Commercial' ! Correct values ! ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23,25,26, & ! 27,31,32,33 /) ! How ever, USGS by default has 24 so... ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23 /) lic_wk = (/ 24 /) oce_wk = (/ 16 /) sic_wk = (/ 24 /) CASE ('MODIFIED_IGBP_MODIS_NOAH') Ncattot = 33 Nter_wk = 28 Nlic_wk = 1 Noce_wk = 1 Nsic_wk = 1 IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk) ALLOCATE(ter_wk(Nter_wk)) IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk) ALLOCATE(lic_wk(Nlic_wk)) IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk) ALLOCATE(oce_wk(Noce_wk)) IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk) ALLOCATE(sic_wk(Nsic_wk)) ! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3) !! ! 1,'Evergreen Needleleaf Forest' 2,'Evergreen Broadleaf Forest' ! 3,'Deciduous Needleleaf Forest' 4,'Deciduous Broadleaf Forest' ! 5,'Mixed Forests' 6,'Closed Shrublands' ! 7,'Open Shrublands' 8,'Woody Savannas' ! 9,'Savannas' 10,'Grasslands' ! 11,'Permanent wetlands' 12,'Croplands' ! 13,'Urban and Built-Up' 14,'cropland/natural vegetation mosaic' ! 15,'Snow and Ice' 16,'Barren or Sparsely Vegetated' ! 17,'Water' 18,'Wooded Tundra' ! 19,'Mixed Tundra' 20,'Barren Tundra' ! 21,'Unassigned' 22,'Unassigned' ! 23,'Unassigned' 24,'Unassigned' ! 25,'Unassigned' 26,'Unassigned' ! 27,'Unassigned' 28,'Unassigned' ! 29,'Unassigned' 30,'Unassigned' ! 31,'Low Intensity Residential ' 32,'High Intensity Residential' ! 33,'Industrial or Commercial' ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14,16,18,19,20,31,32,33 /) lic_wk = (/ 15 /) oce_wk = (/ 17 /) sic_wk = (/ 15 /) CASE ('SiB') Ncattot = 16 Nter_wk = 14 Nlic_wk = 1 Noce_wk = 1 Nsic_wk = 1 IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk) ALLOCATE(ter_wk(Nter_wk)) IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk) ALLOCATE(lic_wk(Nlic_wk)) IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk) ALLOCATE(oce_wk(Noce_wk)) IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk) ALLOCATE(sic_wk(Nsic_wk)) ! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3) !! ! 1,'Evergreen Broadleaf Trees' 2,'Broadleaf Deciduous Trees' ! 3,'Deciduous and Evergreen Trees' 4,'Evergreen Needleleaf Trees' ! 5,'Deciduous Needleleaf Trees' 6,'Ground Cover with Trees and Shrubs' ! 7,'Groundcover Only' 8,'Broadleaf Shrubs with Perennial Ground Cover' ! 9,'Broadleaf Shrubs with Bare Soil' 10,'Groundcover with Dwarf Trees and Shrubs' ! 11,'Bare Soil' 12,'Agriculture or C3 Grassland' ! 13,'Persistent Wetland' 14,'Dry Coastal Complexes' ! 15,'Water' 16,'Ice Cap and Glacier' ter_wk = (/ 1,2,3,4,5,6,7,8,9,10,11,12,13,14 /) lic_wk = (/ 16 /) oce_wk = (/ 15 /) sic_wk = (/ 16 /) CASE ('LW12') Ncattot = 3 Nter_wk = 1 Nlic_wk = 1 Noce_wk = 1 Nsic_wk = 1 IF (ALLOCATED(ter_wk)) DEALLOCATE(ter_wk) ALLOCATE(ter_wk(Nter_wk)) IF (ALLOCATED(lic_wk)) DEALLOCATE(lic_wk) ALLOCATE(lic_wk(Nlic_wk)) IF (ALLOCATED(oce_wk)) DEALLOCATE(oce_wk) ALLOCATE(oce_wk(Noce_wk)) IF (ALLOCATED(sic_wk)) DEALLOCATE(sic_wk) ALLOCATE(sic_wk(Nsic_wk)) ! It could be done much more 'elegantly' ... reading values from LANDUSE.TBL (v3.3) !! ! 1,'Land' 2,'Water' ! 3,'Snow or Ice' ter_wk = (/ 1 /) lic_wk = (/ 3 /) oce_wk = (/ 2 /) sic_wk = (/ 3 /) CASE DEFAULT PRINT *,TRIM(errmsg) PRINT *,' ' // TRIM(fname) // ": landuse data-base '" // TRIM(lndcn) // & "' not ready !!!" STOP END SELECT ! Checking the correct number of categories IF (Ncattot /= ncat) THEN PRINT *,TRIM(wrnmsg) PRINT *,' ' // TRIM(fname) // ": '" // TRIM(lndcn) // "' landuse data-base " & // ' has ',Ncattot,' land-uses types ',ncat,' passed !!!!' ! STOP END IF ! Mimic WRF landuse with 4 LMDZ Ksoils !! IF (ALLOCATED(wk)) DEALLOCATE(wk) ! ALLOCATE(wk(4,ncattot)) ALLOCATE(wk(4,ncat)) wk = .FALSE. DO i=1,Nter_wk wk(1,ter_wk(i)) = .TRUE. END DO DO i=1,Nlic_wk wk(2,lic_wk(i)) = .TRUE. END DO DO i=1,Noce_wk wk(3,oce_wk(i)) = .TRUE. END DO DO i=1,Nsic_wk wk(4,sic_wk(i)) = .TRUE. END DO klks_equal = 0 DO i=1,dx DO j=1,dy ! Using double precission numbers to avoid truncation errors... Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(1,:)) kt(i,j) = Ksoilval(1,1) Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(2,:)) kl(i,j) = Ksoilval(1,1) Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(3,:)) ko(i,j) = Ksoilval(1,1) Ksoilval(1,1) = SUM(lndu(i,:,j),MASK=wk(4,:)) ks(i,j) = Ksoilval(1,1) IF (kl(i,j) == ks(i,j)) klks_equal = klks_equal+1 END DO END DO ! Assigning lic/sic to all that ice/snow according to the percentage of ter/oce ! Is there any other better way? Satellites, do not know what is below the ice ?! IF (klks_equal == dx*dy) THEN ! Assuming that kl=ks=ice DO i=1,dx DO j=1,dy IF ((kt(i,j)+ko(i,j)) /= 0.) THEN Ksoilval(1,1) = kl(i,j)*kt(i,j)/(kt(i,j)+ko(i,j)) kl(i,j) = Ksoilval(1,1) Ksoilval(1,1) = ks(i,j)*ko(i,j)/(kt(i,j)+ko(i,j)) ks(i,j) = Ksoilval(1,1) ELSE Ksoilval(1,1) = kl(i,j)/(kl(i,j)+ks(i,j)) Ksoilval(1,2) = ks(i,j)/(kl(i,j)+ks(i,j)) kl(i,j) = Ksoilval(1,1) ks(i,j) = Ksoilval(1,2) END IF ! Introducing WRF land-sea mask IF (mask(i,j) == 1.) THEN kl(i,j)=kl(i,j)+ks(i,j) ks(i,j)=0. ELSE ks(i,j)=kl(i,j)+ks(i,j) kl(i,j)=0. END IF END DO END DO ELSE PRINT *,' ' // TRIM(fname) // ': land-use data provided different lic & sic !' END IF ! Checking if at all the grid points kt+kl+ko+ks == 1. !! totlndu = kt+kl+ko+ks Nnones=0 DO i=1,dx DO j=1,dy IF (totlndu(i,j) /= 1. ) Nnones = Nnones+1 IF (kt(i,j) + ko(i,j) == 0.)THEN PRINT *,TRIM(wrnmsg) PRINT *,' ' // TRIM(fname) // ': point without land or water category !!' PRINT *,'i j WRF_ter ter ______________' DO k=1,Nter_wk PRINT *,i,j,lndu(i,ter_wk(k),j), kt(i,j) END DO PRINT *,'i j WRF_oce oce ______________' DO k=1,Noce_wk PRINT *,i,j,lndu(i,oce_wk(k),j), ko(i,j) END DO PRINT *,' wrf%landusef: ', lndu(i,:,j) ! STOP END IF END DO END DO ij=1 IF (Nnones /= 0) THEN PRINT *,TRIM(errmsg) PRINT *,' ' // TRIM(fname) // ': ', Nnones,' grid points do not have a ' // & 'TOTlanduse == 1.!!!!!' PRINT *,'pt i j ter+lic+oce+sic : ter lic oce sic ______________' DO i=1,dx DO j=1,dy IF (totlndu(i,j) /= 1.) THEN PRINT *,ij,'; ',i,j,totlndu(i,j),' : ',kt(i,j),kl(i,j),ko(i,j),ks(i,j) ij=ij+1 END IF END DO END DO END IF DEALLOCATE(ter_wk) DEALLOCATE(lic_wk) DEALLOCATE(oce_wk) DEALLOCATE(sic_wk) DEALLOCATE(wk) RETURN END SUBROUTINE WRFlanduse_LMDZKsoil SUBROUTINE interpolate_layers(NvalsA,posA,valsA,NvalsB,posB,valsB) ! Subroutine to interpolate values between different layers (constant values along the depth of a layer) IMPLICIT NONE INTEGER, INTENT(IN) :: NvalsA,NvalsB REAL, DIMENSION(NvalsA), INTENT(IN) :: posA,valsA REAL, DIMENSION(NvalsB), INTENT(IN) :: posB REAL, DIMENSION(NvalsB), INTENT(OUT) :: valsB ! Local INTEGER :: i,k REAL :: xA,yA,xB,yB !!!!!!! Variables ! A: original values used to interpolate ! B: values to obtain ! Nvals[A/B]: number of values ! pos[A/B]: depths of each layer ! vals[A/B]: values of each layer !!!!!!! Functions/Subroutines ! interpolate_lin: subroutine to linearly interpolate i=1 DO k=1,NvalsB IF (posB(k) > posA(i)) i = i+1 IF (i-1 < 1) THEN xA = posA(i) yA = valsA(i) xB = posA(i+1) yB = valsA(i+1) ELSEIF ( i > NvalsA) THEN xA = posA(i-2) yA = valsA(i-2) xB = posA(i-1) yB = valsA(i-1) ELSE xA = posA(i-1) yA = valsA(i-1) xB = posA(i) yB = valsA(i) END IF IF (posB(k) <= xB .AND. i-1 > 1) THEN CALL interpolate_lin(xA,yA,xB,yB,posB(k),valsB(k)) ELSEIF (posB(k) <= xA .AND. i-1 < 1) THEN valsB(k) = yA + (yB - yA)*(posB(k) - xA)/(xB - xA) ELSE valsB(k) = yB + (yB - yA)*(posB(k) - xB)/(xB - xA) END IF ! PRINT *,'i ',i,' Interpolating for: ',posB(k),' with xA: ', xA,' yA: ',yA,' & ', & ! ' xB: ', xB,' yB: ',yB,' --> ',valsB(k) END DO END SUBROUTINE interpolate_layers SUBROUTINE interpolate_lin(x1,y1,x2,y2,x0,y) ! Program to interpolate values y=a+b*x with: (from clWRF modifications) ! a=y1 ! b=(y2-y1)/(x2-x1) ! x=(x2-x0) IMPLICIT NONE REAL, INTENT (IN) :: x1,x2,x0 ! REAL(r8), INTENT (IN) :: y1,y2 ! REAL(r8), INTENT (OUT) :: y ! REAL(r8) :: a,b,x REAL, INTENT (IN) :: y1,y2 REAL, INTENT (OUT) :: y REAL :: a,b,x a=y1 b=(y2-y1)/(x2-x1) IF (x0 .GE. x1) THEN x=x0-x1 ELSE x=x1-x0 b=-b ENDIF y=a+b*x ! PRINT *,'y=a+b*x' ! IF (b .LT. 0.) THEN ! PRINT *,'y=',y1,' ',b,'*x' ! ELSE ! PRINT *,'y=',y1,'+',b,'*x' ! ENDIF ! PRINT *,'y(',x,')=',y END SUBROUTINE interpolate_lin SUBROUTINE table_characteristics(datasetn,tablen,datasetline,Nvar,Nr,Nc) ! Subroutine to read values from a WRF .TBL IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: datasetn CHARACTER(LEN=100), INTENT(IN) :: tablen INTEGER, INTENT(OUT) :: datasetline, Nvar, Nr, Nc ! Local INTEGER :: il, iline, ic CHARACTER(LEN=50) :: fname, errmsg CHARACTER(LEN=200) :: lineval LOGICAL :: existsfile, is_used INTEGER :: funit, Llineval, icoma !!!!!! Variables ! datasetn: name of the data set to use ! tablen: name of the table to read ! Nvar: number of variations (g.e. seasons) ! Nr: table row number ! Nc: table column number ! datasetline: line in the table where the dataset starts fname='table_characteristics' errmsg='ERROR -- error -- ERROR -- error' INQUIRE(FILE=TRIM(tablen), EXIST=existsfile) IF (existsfile) THEN DO funit=10,100 INQUIRE(UNIT=funit, OPENED=is_used) IF (.NOT. is_used) EXIT END DO ELSE PRINT *,TRIM(errmsg) PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!' STOP END IF OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD') ! Looking for the place where data-set starts !! datasetline = -1 iline=1 DO il=1,100000 READ(funit,*,END=100)lineval IF (TRIM(lineval) == TRIM(datasetn)) datasetline = il END DO 100 CONTINUE IF (datasetline == -1) THEN PRINT *,TRIM(errmsg) PRINT *,TRIM(fname),': in file "'// TRIM(tablen) // '" dataset: "'// & TRIM(datasetn) // '" does not exist !!!' STOP END IF ! Looking for number of values !! REWIND(UNIT=funit) DO ic=1,datasetline READ(funit,*)lineval END DO READ(funit,*)Nr,Nvar IF (Nvar > 1) READ(funit,*)lineval READ(funit,'(200A)')lineval Llineval=LEN(TRIM(lineval)) ! PRINT *,' lineval: ******' // TRIM(lineval) // '******', Llineval Nc = 0 ic = 1 DO WHILE (ic <= Llineval) ! PRINT *,lineval(ic:Llineval) icoma = SCAN(lineval(ic:Llineval),',') IF (icoma > 1) THEN Nc = Nc + 1 ! PRINT *,' Nc: ',Nc,' ic ',ic,' icoma: ',icoma ic = ic + icoma ELSE ic = ic + Llineval END IF END DO PRINT *,' dataset :"' // TRIM(datasetn) // '" at line ',datasetline, & ' in table :"'// TRIM(tablen) // '" has N variations: ',Nvar, & ' Number of types: ',Nr, ' Number of values: ',Nc CLOSE(funit) RETURN END SUBROUTINE table_characteristics SUBROUTINE read_table(datasetn,tablen,datasetline,Nvar,Nr,Nc,vals,valns) ! Subroutine to read values from a WRF .TBL IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: datasetn CHARACTER(LEN=100), INTENT(IN) :: tablen INTEGER, INTENT(IN) :: datasetline, Nvar, Nr, Nc REAL, DIMENSION(Nvar,Nr,Nc), INTENT(OUT) :: vals CHARACTER(LEN=50), DIMENSION(Nr), INTENT(OUT) :: valns ! Local INTEGER :: iv, ir, ic CHARACTER(LEN=50) :: fname, errmsg, line LOGICAL :: existsfile, is_used INTEGER :: funit !!!!!!! Variables ! datasetn: name of the data set to use ! tablen: name of the table to read ! Nvar: number of variations (g.e. seasons) ! Nr: table row number ! Nc: table column number ! datasetline: line in the table where the dataset starts ! vals: values of the data-set ! valns: assigned name of each value in the dataset fname='read_table' errmsg='ERROR -- error -- ERROR -- error' INQUIRE(FILE=TRIM(tablen), EXIST=existsfile) IF (existsfile) THEN DO funit=10,100 INQUIRE(UNIT=funit, OPENED=is_used) IF (.NOT. is_used) EXIT END DO ELSE PRINT *,TRIM(errmsg) PRINT *,TRIM(fname),': file "'// TRIM(tablen) // '" does not exist !!!' STOP END IF OPEN(UNIT=funit, FILE=TRIM(tablen), STATUS='OLD') DO ic=1,datasetline+1 READ(funit,*)line END DO DO iv=1, Nvar IF (Nvar > 1) READ(funit,*)line DO ir=1, Nr READ(funit,*)(vals(iv,ir,ic), ic=1,Nc),valns(ir) END DO END DO ! PRINT *,' values dataset :"' // TRIM(datasetn) // '" in table :"' //TRIM(tablen) ! DO iv=1,Nvar ! PRINT *,' variation: ',iv ! DO ir=1,Nr ! PRINT *,' ',(vals(iv,ir,ic), ic=1,Nc),valns(ir) ! END DO ! END DO CLOSE(funit) RETURN END SUBROUTINE read_table SUBROUTINE compute_soil_mass(Nvar,Nr,Nc,vals,valns,soilt,Nnost,NOsoilt,smois,ddx, & ddy,ddz,soilm,soilq) ! Subroutine to compute the total soil mass (kg) from a point with different ! percentages of soil types IMPLICIT NONE INTEGER, INTENT(IN) :: Nvar, Nr, Nc, Nnost REAL, DIMENSION(Nvar,Nr,Nc), INTENT(IN) :: vals CHARACTER(LEN=50), DIMENSION(Nr), INTENT(IN) :: valns REAL, DIMENSION(Nr), INTENT(IN) :: soilt CHARACTER(LEN=50), DIMENSION(Nnost), INTENT(IN) :: NOsoilt REAL, INTENT(IN) :: smois, ddx, ddy, ddz REAL, INTENT(OUT) :: soilm,soilq ! Local INTEGER :: iv, ir, ic CHARACTER(LEN=50) :: fname, errmsg LOGICAL :: is_land !!!!!!! Variables ! Nvar: number of variations (g.e. seasons) ! Nr: table row number ! Nc: table column number ! vals: values of the data-set (col2: density [g cm-3]) ! valns: assigned name of each value in the dataset ! soilt: vector with the percentages of soil types ! Nnost: number of soil types which are not land ! NOsoilt: labels of the non-land soil types ! smois: soil moisture [m3 m-3] ! ddx,ddy,ddz: volume of the grid cell [m] ! soilm: total soil mass [kg] ! soilq: total soil water [kg kg-1] fname='read_table' errmsg='ERROR -- error -- ERROR -- error' soilm = 0. soilq = 0. ! PRINT *,' soilmoisture: ',smois DO ir=1,Nr is_land = .TRUE. IF (soilt(ir) /= 0.) THEN DO iv=1, Nnost IF (TRIM(valns(ir)) == TRIM(NOsoilt(iv))) is_land = .FALSE. END DO IF (is_land) THEN ! PRINT *,' soil moisture range: ', vals(1,ir,3), vals(1,ir,5) soilm = soilm + soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000. soilq = soilq + soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000. ! PRINT *,' Percentage of soil "' // TRIM(valns(ir)) //'": ',soilt(ir)*100.,& ! ' % density: ',vals(1,ir,2),' partial mass: ', & ! soilt(ir)*vals(1,ir,2)*ddx*ddy*ddz*1000. , ' partial water: ', & ! soilt(ir)*vals(1,ir,2)*smois*ddx*ddy*ddz*1000. END IF END IF END DO ! Remember 1 kg H20 = 1 l = 1mm / m2 !! PRINT *,' TOTAL. soil mass: ',soilm,' [kg] soil water: ',soilq,' [kg]', & !! ' equiv. water depth in soil ',soilq/(ddx*ddy*1000.),' [m]' END SUBROUTINE compute_soil_mass SUBROUTINE compute_TOTsoil_mass_water(datname, tabname, Nsoiltypes, soiltypest, & soiltypesb, Nsoillayers, smois, soillayers, sizeX, sizeY, soilTOTmass, & soilTOThumidity) ! Subroutine to compute total soil mass and water from a given grid point using a ! given data-set from SOILPARM.TBL IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: datname CHARACTER(LEN=100), INTENT(IN) :: tabname INTEGER, INTENT(IN) :: Nsoiltypes, Nsoillayers REAL, DIMENSION(Nsoiltypes), INTENT(IN) :: soiltypest, soiltypesb REAL, DIMENSION(Nsoillayers), INTENT(IN) :: smois, soillayers REAL, INTENT(IN) :: sizeX, sizeY REAL, INTENT(OUT) :: soilTOTmass, & soilTOThumidity ! Local INTEGER :: ix,iy,iz REAL, ALLOCATABLE, DIMENSION(:,:,:) :: values CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:) :: namevalues,NOsoiltypes INTEGER :: linedataset,Nvariations, & Nrows, Ncols INTEGER :: NNOsoiltypes REAL, ALLOCATABLE, DIMENSION(:) :: soiltypes REAL :: dz,soilmass,soilhumidity CHARACTER(LEN=50) :: fname, errmsg !!!!!!! Variables ! datname: name of the data set to use ! tabname: name of the table to read ! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons), ! table row (Nrows), table column (Ncols) ! linedataset: line inside data-set where values start ! namevalues: string assigned to each table row ! Nsoiltypes: number of soil types ! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom ! soiltypes: percentage of the soil types at a given grid point (lineray interpolated) ! NOsoiltypes: vector with the names of soil types which are not land fname='"compute_TOTsoil_mass_water"' errmsg='ERROR -- error -- ERROR -- error' CALL table_characteristics(datname, tabname, linedataset, Nvariations, & Nrows, Ncols) IF (Nrows /= Nsoiltypes) THEN PRINT *,errmsg PRINT ' ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes, & ' does not concide with the number of values: ', Nrows, & ' from the data-set "' // TRIM(datname) // '" !!!' STOP END IF IF (ALLOCATED(values)) DEALLOCATE(values) ALLOCATE(values(Nvariations,Nrows,Ncols)) IF (ALLOCATED(namevalues)) DEALLOCATE(namevalues) ALLOCATE(namevalues(Nrows)) CALL read_table(datname, tabname, linedataset, Nvariations, Nrows, Ncols, & values, namevalues) IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes) ALLOCATE(soiltypes(Nrows)) ! On WRF3.3 SOILPARM.TBL SELECT CASE (TRIM(datname)) CASE ('STAS') NNOsoiltypes = 1 IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes) ALLOCATE(NOsoiltypes(NNOsoiltypes)) NOsoiltypes=(/ 'WATER' /) CASE ('STAS-RUC') NNOsoiltypes = 1 IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes) ALLOCATE(NOsoiltypes( NNOsoiltypes)) NOsoiltypes=(/ 'WATER' /) END SELECT soilTOTmass = 0. soilTOThumidity = 0. DO iz=1,Nsoillayers IF (iz == 1 ) THEN dz = soillayers(iz) ELSE dz = soillayers(iz) - soillayers(iz - 1) END IF IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN DO ix=1,Nrows CALL interpolate_lin(soillayers(1), soiltypest(ix), & soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix)) END DO ELSEIF (iz == 1) THEN soiltypes = soiltypest ELSE soiltypes = soiltypesb END IF ! PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes, & NNOsoiltypes,NOsoiltypes,smois(iz),1.,1.,dz,soilmass,soilhumidity) soilTOTmass = soilTOTmass + soilmass soilTOThumidity = soilTOThumidity + soilhumidity END DO DEALLOCATE(values) DEALLOCATE(namevalues) DEALLOCATE(soiltypes) DEALLOCATE(NOsoiltypes) END SUBROUTINE compute_TOTsoil_mass_water SUBROUTINE compute_TOTsoil_mass_water_values(datname, Nvariations, Nrows, Ncols, & values, namevalues, ip, jp, Nsoiltypes, soiltypest, soiltypesb, Nsoillayers, & smois, soillayers, sizeX, sizeY, soilTOTmass, soilTOThumidity) ! Subroutine to compute total soil mass and water from a given grid point using the ! values from a given data-set of SOILPARM.TBL IMPLICIT NONE CHARACTER(LEN=50), INTENT(IN) :: datname INTEGER, INTENT(IN) :: Nvariations, Nrows, Ncols REAL, DIMENSION(Nvariations, Nrows, Ncols), INTENT(IN) :: values CHARACTER(LEN=50), DIMENSION(Nrows), INTENT(IN) :: namevalues INTEGER, INTENT(IN) :: ip, jp INTEGER, INTENT(IN) :: Nsoiltypes, Nsoillayers REAL, DIMENSION(Nsoiltypes), INTENT(IN) :: soiltypest, soiltypesb REAL, DIMENSION(Nsoillayers), INTENT(IN) :: smois, soillayers REAL, INTENT(IN) :: sizeX, sizeY REAL, INTENT(OUT) :: soilTOTmass, & soilTOThumidity ! Local INTEGER :: ix,iy,iz CHARACTER(LEN=50), ALLOCATABLE, DIMENSION(:) :: NOsoiltypes INTEGER :: linedataset INTEGER :: NNOsoiltypes REAL, ALLOCATABLE, DIMENSION(:) :: soiltypes REAL :: dz,soilmass,soilhumidity CHARACTER(LEN=50) :: fname, errmsg !!!!!!! Variables ! datname: name of the data set to use ! values: mmatrix with the name of the values: variation (Nvariations, g.e. seasons), ! table row (Nrows), table column (Ncols) ! linedataset: line inside data-set where values start ! namevalues: string assigned to each table row ! ip, jp: coordinates of the grid point ! Nsoiltypes: number of soil types ! soiltypes[t/b]: percentage of the soil types at a given grid point on top/bottom ! soiltypes: percentage of the soil types at a given grid point (lineray interpolated) ! NOsoiltypes: vector with the names of soil types which are not land fname='"compute_TOTsoil_mass_water_values"' errmsg='ERROR -- error -- ERROR -- error' IF (Nrows /= Nsoiltypes) THEN PRINT *,errmsg PRINT ' ' // TRIM(fname) // ': Number of soil-types: ',Nsoiltypes, & ' does not concide with the number of values: ', Nrows, & ' from the data-set "' // TRIM(datname) // '" !!!' STOP END IF IF (ALLOCATED(soiltypes)) DEALLOCATE(soiltypes) ALLOCATE(soiltypes(Nrows)) ! On WRF3.3 SOILPARM.TBL SELECT CASE (TRIM(datname)) CASE ('STAS') NNOsoiltypes = 1 IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes) ALLOCATE(NOsoiltypes(NNOsoiltypes)) NOsoiltypes=(/ 'WATER' /) CASE ('STAS-RUC') NNOsoiltypes = 1 IF (ALLOCATED(NOsoiltypes)) DEALLOCATE(NOsoiltypes) ALLOCATE(NOsoiltypes( NNOsoiltypes)) NOsoiltypes=(/ 'WATER' /) END SELECT soilTOTmass = 0. soilTOThumidity = 0. DO iz=1,Nsoillayers IF (iz == 1 ) THEN dz = soillayers(iz) ELSE dz = soillayers(iz) - soillayers(iz - 1) END IF IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN DO ix=1,Nrows CALL interpolate_lin(soillayers(1), soiltypest(ix), & soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz), soiltypes(ix)) END DO ELSEIF (iz == 1) THEN soiltypes = soiltypest ELSE soiltypes = soiltypesb END IF ! PRINT *,'lev: ',iz,' dz: ',dz,' soil types: ' ,soiltypes CALL compute_soil_mass(Nvariations,Nrows,Ncols,values,namevalues,soiltypes, & NNOsoiltypes,NOsoiltypes,smois(iz),sizeX,sizeY,dz,soilmass,soilhumidity) soilTOTmass = soilTOTmass + soilmass soilTOThumidity = soilTOThumidity + soilhumidity END DO IF ( (soilTOTmass < 0.) .OR. (soilTOThumidity < 0.) ) THEN PRINT *,errmsg PRINT *,' ' // TRIM(fname) // ' at point ',ip, ', ', jp, & ' negative total soil mass: ',soilTOTmass, ' or soil water content', & soilTOThumidity PRINT *,' Soil_name density soil_types_________' DO ix=1, Nrows PRINT *,ix, TRIM(namevalues(ix)), values(1,ix,2), soiltypes(ix) END DO PRINT *,' depth_layer partial_mass soil+moisure partial_moisture__________' DO iz=1,Nsoillayers IF (iz == 1 ) THEN dz = soillayers(iz) ELSE dz = soillayers(iz) - soillayers(iz - 1) END IF IF (( iz > 1) .AND. (iz < Nsoillayers)) THEN DO ix=1,Nrows CALL interpolate_lin(soillayers(1), soiltypest(ix), & soillayers(Nsoillayers), soiltypesb(ix), soillayers(iz),soiltypes(ix)) END DO ELSEIF (iz == 1) THEN soiltypes = soiltypest ELSE soiltypes = soiltypesb END IF PRINT *,' ',iz, dz, values(1,iz,2)*sizeX*sizeY*dz*1000., smois(iz), & smois(iz)*values(1,iz,2)*sizeX*sizeY*dz*1000. END DO STOP END IF DEALLOCATE(soiltypes) DEALLOCATE(NOsoiltypes) END SUBROUTINE compute_TOTsoil_mass_water_values ! cat Registry_out.LMDZ | awk '{print " pwrf_grid%"tolower($3)" "$4" = p"tolower($3)" "$4}' >> get_lmdz_out.f90 ! cat get_lmdz_out.f90n' ',' '{print $3}' | tr '\n SUBROUTINE get_lmdz_out(ddx, ddy, ddz, bdy, dlmdz, Nks, Nz0, nqtot, ivap, iliq, & iflag_pbl, iflag_con, iflag_thermals, iflag_wake, read_climoz, itap, & pdtphys, paprs, pplay, pphi, pphis, presnivs, qx, & d_u, d_v, d_t, d_qx, d_ps, & pwrf_grid) ! Subroutine to get all the LMDZ output variables into the WRF Registry structure ! WRF modules USE module_domain_type ! LMDZ modules USE indice_sol_mod USE phys_state_var_mod USE phys_output_var_mod ! USE pbl_surface_mod USE phys_local_var_mod USE comgeomphy ! WRF_LMDZ modules USE output_lmdz_NOmodule ! netcdf USE netcdf, ONLY: nf90_fill_real IMPLICIT NONE TYPE(domain) , TARGET :: pwrf_grid INTEGER, INTENT(IN) :: ddx, ddy, ddz, bdy, & dlmdz, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals, & iflag_wake, read_climoz, itap REAL, INTENT(IN) :: pdtphys REAL, DIMENSION(dlmdz), INTENT(IN) :: pphis, d_ps REAL, DIMENSION(ddz), INTENT(IN) :: presnivs REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN) :: paprs REAL, DIMENSION(dlmdz,ddz), INTENT(IN) :: pplay, pphi, d_u, d_v, d_t REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN) :: qx, d_qx ! EXTERNAL suphel ! Local INTEGER :: dx, dy, ix, iy, ixy, iz INTEGER :: i,j,k,l,n INTEGER :: nsrf REAL, DIMENSION(dlmdz,ddz) :: var3d REAL, DIMENSION(dlmdz,ddz+1) :: var3dstagg REAL :: RDAY, RG, RD, RMO3 CHARACTER(LEN=50) :: fname, errmsg CHARACTER(LEN=4) :: bb2 REAL, DIMENSION(dlmdz) :: zx_tmp_fi2d ! variable temporaire grille physique REAL, DIMENSION(dlmdz,ddz) :: zx_tmp_fi3d ! variable temporaire pour champs 3D REAL, DIMENSION(dlmdz,ddz+1) :: zx_tmp_fi3d1 !variable temporaire pour champs 3D (kelvp1) REAL(KIND=8), DIMENSION(dlmdz,ddz) :: zx_tmp2_fi3d ! variable temporaire pour champs 3D REAL, PARAMETER :: dobson_u = 2.1415e-05 ! klon=dldmz ! klev=ddz !! To think about / discuss with Frederic INCLUDE "declare_STDlev.h" ! INCLUDE "calcul_STDlev.h" !!!!!!! Variables ! pwrf_grid: WRF grid structure ! L. Fita, LMD January 2014. It should work using suphel, but not? ! CALL suphel fname='"get_lmdz_out"' errmsg='ERROR -- error -- ERROR -- error' RDAY = 86400. RG = 9.81 RD = 287. RMO3=47.9942 ! Remove variable assignation already done in get_lmdz? These variables are not used ! in restart/boundaries DO ix=1,ddx DO iy=1,ddy ixy=ddx*(iy-1)+ix pwrf_grid%lages_ter(ix,iy) = agesno(ixy,1) pwrf_grid%lages_lic(ix,iy) = agesno(ixy,2) pwrf_grid%lages_oce(ix,iy) = agesno(ixy,3) pwrf_grid%lages_sic(ix,iy) = agesno(ixy,4) !! pwrf_grid%lahyb(misc??) = pahyb misc pwrf_grid%laire(ix,iy) = airephy(ixy) pwrf_grid%laireter(ix,iy) = paire_ter(ixy) pwrf_grid%lalb1(ix,iy) = albsol1(ixy) pwrf_grid%lalb2(ix,iy) = albsol2(ixy) pwrf_grid%lalbe_ter(ix,iy) = falb1(ixy,1) pwrf_grid%lalbe_lic(ix,iy) = falb1(ixy,2) pwrf_grid%lalbe_oce(ix,iy) = falb1(ixy,3) pwrf_grid%lalbe_sic(ix,iy) = falb1(ixy,4) pwrf_grid%lale(ix,iy) = ale(ixy) pwrf_grid%lale_bl(ix,iy) = ale_bl(ixy) pwrf_grid%lale_wk(ix,iy) = ale_wake(ixy) pwrf_grid%lalp(ix,iy) = alp(ixy) pwrf_grid%lalp_bl(ix,iy) = alp_bl(ixy) pwrf_grid%lalp_wk(ix,iy) = alp_wake(ixy) !! pwrf_grid%lalt(misc??) = palt misc !! pwrf_grid%lbhyb(misc??) = pbhyb misc pwrf_grid%lbils(ix,iy) = bils(ixy) pwrf_grid%lbils_diss(ix,iy) = bils_diss(ixy) pwrf_grid%lbils_ec(ix,iy) = bils_ec(ixy) pwrf_grid%lbils_enthalp(ix,iy) = bils_enthalp(ixy) pwrf_grid%lbils_kinetic(ix,iy) = bils_kinetic(ixy) pwrf_grid%lbils_latent(ix,iy) = bils_latent(ixy) pwrf_grid%lbils_tke(ix,iy) = bils_tke(ixy) pwrf_grid%lcape(ix,iy) = cape(ixy) pwrf_grid%lcape_max(ix,iy) = cape(ixy) pwrf_grid%lcdrh(ix,iy) = cdragh(ixy) pwrf_grid%lcdrm(ix,iy) = cdragm(ixy) pwrf_grid%lcldh(ix,iy) = cldh(ixy) pwrf_grid%lcldl(ix,iy) = cldl(ixy) pwrf_grid%lcldm(ix,iy) = cldm(ixy) pwrf_grid%lcldq(ix,iy) = cldq(ixy) pwrf_grid%lcldt(ix,iy) = cldt(ixy) pwrf_grid%lcontfracatm(ix,iy) = pctsrf(ixy,is_ter)+pctsrf(ixy,is_lic) pwrf_grid%lcontfracor(ix,iy) = pctsrf(ixy,is_ter) !NOTknown! pwrf_grid%lcumpb(ix,iy) = cumpb(ixy) !NOTknown! pwrf_grid%lcumrn(ix,iy) = cumrn(ixy) pwrf_grid%ldthmin(ix,iy) = dthmin(ixy) pwrf_grid%ldtsvdft(ix,iy) = d_ts(ixy,is_ter) pwrf_grid%ldtsvdfo(ix,iy) = d_ts(ixy,is_oce) pwrf_grid%ldtsvdfg(ix,iy) = d_ts(ixy,is_lic) pwrf_grid%ldtsvdfi(ix,iy) = d_ts(ixy,is_sic) ! L. Fita, October 2014. with fevap there is not evaporation! ! pwrf_grid%levap(ix,iy) = fevap(ixy,1) pwrf_grid%levap(ix,iy) = evap(ixy) pwrf_grid%levap_ter(ix,iy) = fevap(ixy,1) pwrf_grid%levap_lic(ix,iy) = fevap(ixy,2) pwrf_grid%levap_oce(ix,iy) = fevap(ixy,3) pwrf_grid%levap_sic(ix,iy) = fevap(ixy,4) pwrf_grid%levappot_ter(ix,iy) = evap_pot(ixy,1) pwrf_grid%levappot_lic(ix,iy) = evap_pot(ixy,2) pwrf_grid%levappot_oce(ix,iy) = evap_pot(ixy,3) pwrf_grid%levappot_sic(ix,iy) = evap_pot(ixy,4) pwrf_grid%lfbase(ix,iy) = ema_cbmf(ixy) pwrf_grid%lfder(ix,iy) = fder(ixy) pwrf_grid%lffonte(ix,iy) = zxffonte(ixy) pwrf_grid%lflat(ix,iy) = zxfluxlat(ixy) pwrf_grid%lflw_ter(ix,iy) = fsollw(ixy,1) pwrf_grid%lflw_lic(ix,iy) = fsollw(ixy,2) pwrf_grid%lflw_oce(ix,iy) = fsollw(ixy,3) pwrf_grid%lflw_sic(ix,iy) = fsollw(ixy,4) pwrf_grid%lfqcalving(ix,iy) = zxfqcalving(ixy) pwrf_grid%lfqfonte(ix,iy) = zxfqfonte(ixy) pwrf_grid%lfract_ter(ix,iy) = pctsrf(ixy,1) pwrf_grid%lfract_lic(ix,iy) = pctsrf(ixy,2) pwrf_grid%lfract_oce(ix,iy) = pctsrf(ixy,3) pwrf_grid%lfract_sic(ix,iy) = pctsrf(ixy,4) pwrf_grid%lfsnow(ix,iy) = zfra_o(ixy) pwrf_grid%lfsw_ter(ix,iy) = fsolw(ixy,1) pwrf_grid%lfsw_lic(ix,iy) = fsolw(ixy,2) pwrf_grid%lfsw_oce(ix,iy) = fsolw(ixy,3) pwrf_grid%lfsw_sic(ix,iy) = fsolw(ixy,4) !NOtnow! pwrf_grid%lftime_con(ix,iy) = float(itau_con)/float(itap) pwrf_grid%lftime_th(ix,iy) = 0. pwrf_grid%liwp(ix,iy) = fiwp(ixy) !! pwrf_grid%llat j = pllat j pwrf_grid%llat_ter(ix,iy) = fluxlat(ixy,1) pwrf_grid%llat_lic(ix,iy) = fluxlat(ixy,2) pwrf_grid%llat_oce(ix,iy) = fluxlat(ixy,3) pwrf_grid%llat_sic(ix,iy) = fluxlat(ixy,4) pwrf_grid%llmaxth(ix,iy) = lmax_th(ixy) !! pwrf_grid%llon i = pllon i pwrf_grid%llwdn200(ix,iy) = lwdn200(ixy) pwrf_grid%llwdn200clr(ix,iy) = lwdn200clr(ixy) pwrf_grid%llwdnsfc(ix,iy) = sollwdown(ixy) pwrf_grid%llwdnsfcclr(ix,iy) = sollwdownclr(ixy) pwrf_grid%llwdownor(ix,iy) = sollwdown(ixy) pwrf_grid%llwp(ix,iy) = flwp(ixy) pwrf_grid%llwup200(ix,iy) = lwup200(ixy) pwrf_grid%llwup200clr(ix,iy) = lwup200clr(ixy) pwrf_grid%llwupsfc(ix,iy) = sollwdown(ixy)-sollw(ixy) pwrf_grid%llwupsfcclr(ix,iy) = -1.*lwdn0(ixy,1)-sollw0(ixy) pwrf_grid%lmsnow(ix,iy) = snow_o(ixy) pwrf_grid%lndayrain(ix,iy) = nday_rain(ixy) pwrf_grid%lnettop(ix,iy) = topsw(ixy)-toplw(ixy) pwrf_grid%lpbase(ix,iy) = ema_pcb(ixy) pwrf_grid%lphis(ix,iy) = pphis(ixy) pwrf_grid%lplcl(ix,iy) = plcl(ixy) pwrf_grid%lplfc(ix,iy) = plfc(ixy) pwrf_grid%lpluc(ix,iy) = rain_con(ixy) + snow_con(ixy) pwrf_grid%lplul(ix,iy) = rain_lsc(ixy) + snow_lsc(ixy) pwrf_grid%lplulst(ix,iy) = plul_st(ixy) pwrf_grid%lplulth(ix,iy) = plul_th(ixy) pwrf_grid%lpourc_ter(ix,iy) = pctsrf(ixy,1)*100. pwrf_grid%lpourc_lic(ix,iy) = pctsrf(ixy,2)*100. pwrf_grid%lpourc_oce(ix,iy) = pctsrf(ixy,3)*100. pwrf_grid%lpourc_sic(ix,iy) = pctsrf(ixy,4)*100. pwrf_grid%lprecip(ix,iy) = rain_fall(ixy) + snow_fall(ixy) !! pwrf_grid%lpresnivs k = plpresnivs k pwrf_grid%lprw(ix,iy) = prw(ixy) pwrf_grid%lpsol(ix,iy) = paprs(ixy,1) pwrf_grid%lptop(ix,iy) = ema_pct(ixy) pwrf_grid%lq2m(ix,iy) = zq2m(ixy) pwrf_grid%lqsat2m(ix,iy) = qsat2m(ixy) pwrf_grid%lqsol(ix,iy) = qsol(ixy) pwrf_grid%lqsurf(ix,iy) = zxqsurf(ixy) pwrf_grid%lradsol(ix,iy) = radsol(ixy) pwrf_grid%lrh2m(ix,iy) = MIN(100.,rh2m(ixy)*100.) pwrf_grid%lrh2m_max(ix,iy) = MIN(100.,rh2m(ixy)*100.) pwrf_grid%lrh2m_min(ix,iy) = MIN(100.,rh2m(ixy)*100.) pwrf_grid%lrugs(ix,iy) = zxrugs(ixy) pwrf_grid%lrugs_ter(ix,iy) = frugs(ixy,1) pwrf_grid%lrugs_lic(ix,iy) = frugs(ixy,2) pwrf_grid%lrugs_oce(ix,iy) = frugs(ixy,3) pwrf_grid%lrugs_sic(ix,iy) = frugs(ixy,4) pwrf_grid%lsens(ix,iy) = -1.*sens(ixy) pwrf_grid%lsicf(ix,iy) = pctsrf(ixy,is_sic) pwrf_grid%ls_lcl(ix,iy) = s_lcl(ixy) pwrf_grid%lslp(ix,iy) = slp(ixy) pwrf_grid%lsnow(ix,iy) = snow_fall(ixy) pwrf_grid%lsnowl(ix,iy) = snow_lsc(ixy) pwrf_grid%lsoll(ix,iy) = sollw(ixy) pwrf_grid%lsoll0(ix,iy) = sollw0(ixy) pwrf_grid%lsollwdown(ix,iy) = sollwdown(ixy) pwrf_grid%lsols(ix,iy) = solsw(ixy) pwrf_grid%lsols0(ix,iy) = solsw0(ixy) pwrf_grid%ls_pblh(ix,iy) = s_pblh(ixy) pwrf_grid%ls_pblt(ix,iy) = s_pblt(ixy) pwrf_grid%ls_therm(ix,iy) = s_therm(ixy) pwrf_grid%lswdn200(ix,iy) = swdn200(ixy) pwrf_grid%lswdn200clr(ix,iy) = swdn200clr(ixy) pwrf_grid%lswdnsfc(ix,iy) = swdn(ixy,1) pwrf_grid%lswdnsfcclr(ix,iy) = swdn0(ixy,1) pwrf_grid%lswdntoa(ix,iy) = swdn(ixy,klev+1) pwrf_grid%lswdntoaclr(ix,iy) = swdn0(ixy,klev+1) pwrf_grid%lswdownor(ix,iy) = solsw(ixy)/(1.-albsol1(ixy)) pwrf_grid%lswnetor(ix,iy) = fsolsw(ixy,is_ter) pwrf_grid%lswup200(ix,iy) = swup200(ixy) pwrf_grid%lswup200clr(ix,iy) = swup200clr(ixy) pwrf_grid%lswupsfc(ix,iy) = swup(ixy,1) pwrf_grid%lswupsfcclr(ix,iy) = swup0(ixy,1) pwrf_grid%lswuptoa(ix,iy) = swup(ixy,klev+1) pwrf_grid%lswuptoaclr(ix,iy) = swup0(ixy,klev+1) pwrf_grid%lt2m(ix,iy) = zt2m(ixy) pwrf_grid%lt2m_max(ix,iy) = zt2m(ixy) pwrf_grid%lt2m_min(ix,iy) = zt2m(ixy) pwrf_grid%lt2m_ter(ix,iy) = t2m(ixy,1) pwrf_grid%lt2m_lic(ix,iy) = t2m(ixy,2) pwrf_grid%lt2m_oce(ix,iy) = t2m(ixy,3) pwrf_grid%lt2m_sic(ix,iy) = t2m(ixy,4) pwrf_grid%ltaux(ix,iy) = 0. pwrf_grid%ltauy(ix,iy) = 0. DO iz=1,4 pwrf_grid%ltaux(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)* & fluxu(ixy,1,iz) pwrf_grid%ltauy(ix,iy) = pwrf_grid%ltauy(ix,iy) + pctsrf(ixy,iz)* & fluxv(ixy,1,iz) END DO pwrf_grid%ltaux_ter(ix,iy) = fluxu(ixy,1,1) pwrf_grid%ltaux_lic(ix,iy) = fluxu(ixy,1,2) pwrf_grid%ltaux_oce(ix,iy) = fluxu(ixy,1,3) pwrf_grid%ltaux_sic(ix,iy) = fluxu(ixy,1,4) pwrf_grid%ltauy_ter(ix,iy) = fluxv(ixy,1,1) pwrf_grid%ltauy_lic(ix,iy) = fluxv(ixy,1,2) pwrf_grid%ltauy_oce(ix,iy) = fluxv(ixy,1,3) pwrf_grid%ltauy_sic(ix,iy) = fluxv(ixy,1,4) !! pwrf_grid%ltime - = pltime - !! pwrf_grid%ltime_bounds {lbnds} = pltime_bounds {lbnds} !! pwrf_grid%lti86400 - = plti86400 - IF ( (pctsrf(ixy,is_oce).GT.epsfra).OR.(pctsrf(ixy,is_sic).GT.epsfra) ) THEN pwrf_grid%lt_oce_sic(ix,iy) = (ftsol(ixy, is_oce) * pctsrf(ixy,is_oce)+ & ftsol(ixy, is_sic)*pctsrf(ixy,is_sic))/(pctsrf(ixy,is_oce)+ & pctsrf(ixy,is_sic)) ELSE pwrf_grid%lt_oce_sic(ix,iy) = 273.15 END IF !! pwrf_grid%lt_op_00001800(misc??) = pt_op_00001800 misc !! pwrf_grid%lt_op_00001800_bnds(misc??) = pt_op_00001800_bnds misc pwrf_grid%ltopl(ix,iy) = toplw(ixy) pwrf_grid%ltopl0(ix,iy) = toplw0(ixy) pwrf_grid%ltops(ix,iy) = topsw(ixy) pwrf_grid%ltops0(ix,iy) = topsw0(ixy) pwrf_grid%ltpot(ix,iy) = tpot(ixy) pwrf_grid%ltpote(ix,iy) = tpote(ixy) pwrf_grid%ltsol(ix,iy) = zxtsol(ixy) pwrf_grid%ltsol_ter(ix,iy) = ftsol(ixy,1) pwrf_grid%ltsol_lic(ix,iy) = ftsol(ixy,2) pwrf_grid%ltsol_oce(ix,iy) = ftsol(ixy,3) pwrf_grid%ltsol_sic(ix,iy) = ftsol(ixy,4) pwrf_grid%lu10m(ix,iy) = zu10m(ixy) pwrf_grid%lu10m_ter(ix,iy) = u10m(ixy,1) pwrf_grid%lu10m_lic(ix,iy) = u10m(ixy,2) pwrf_grid%lu10m_oce(ix,iy) = u10m(ixy,3) pwrf_grid%lu10m_sic(ix,iy) = u10m(ixy,4) pwrf_grid%lue(ix,iy) = ue(ixy) pwrf_grid%luq(ix,iy) = uq(ixy) pwrf_grid%lustar(ix,iy) = zustar(ixy) pwrf_grid%lustar_ter(ix,iy) = ustar(ixy,1) pwrf_grid%lustar_lic(ix,iy) = ustar(ixy,2) pwrf_grid%lustar_oce(ix,iy) = ustar(ixy,3) pwrf_grid%lustar_sic(ix,iy) = ustar(ixy,4) pwrf_grid%lv10m(ix,iy) = zv10m(ixy) pwrf_grid%lv10m_ter(ix,iy) = v10m(ixy,1) pwrf_grid%lv10m_lic(ix,iy) = v10m(ixy,2) pwrf_grid%lv10m_oce(ix,iy) = v10m(ixy,3) pwrf_grid%lv10m_sic(ix,iy) = v10m(ixy,4) pwrf_grid%lve(ix,iy) = ve(ixy) pwrf_grid%lvq(ix,iy) = vq(ixy) pwrf_grid%lwake_h(ix,iy) = wake_h(ixy) pwrf_grid%lwake_s(ix,iy) = wake_s(ixy) pwrf_grid%lwape(ix,iy) = wake_pe(ixy) pwrf_grid%lwbeff(ix,iy) = wbeff(ixy) pwrf_grid%lwbilo_ter(ix,iy) = wfbilo(ixy,1) pwrf_grid%lwbilo_lic(ix,iy) = wfbilo(ixy,2) pwrf_grid%lwbilo_oce(ix,iy) = wfbilo(ixy,3) pwrf_grid%lwbilo_sic(ix,iy) = wfbilo(ixy,4) pwrf_grid%lwbils_ter(ix,iy) = wfbils(ixy,1) pwrf_grid%lwbils_lic(ix,iy) = wfbils(ixy,2) pwrf_grid%lwbils_oce(ix,iy) = wfbils(ixy,3) pwrf_grid%lwbils_sic(ix,iy) = wfbils(ixy,4) pwrf_grid%lweakinv(ix,iy) = weak_inversion(ixy) pwrf_grid%lwind10m(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy)) pwrf_grid%lwind10max(ix,iy) = SQRT(zu10m(ixy)*zu10m(ixy)+zv10m(ixy)*zv10m(ixy)) pwrf_grid%lzmax_th(ix,iy) = zmax_th(ixy) ! Output at the standard levels !! ! DO iz=1, nlevSTD ! bb2=clevSTD(iz) ! IF (TRIM(bb2) == "850") THEN ! pwrf_grid%lq850(ix,iy) = q850(ixy) ! pwrf_grid%lt850(ix,iy) = t850(ixy) ! pwrf_grid%lu850(ix,iy) = u850(ixy) ! pwrf_grid%lv850(ix,iy) = v850(ixy) ! pwrf_grid%lw850(ix,iy) = w850(ixy) ! pwrf_grid%lz850(ix,iy) = z850(ixy) ! ELSE IF (TRIM(bb2) == "700") THEN ! pwrf_grid%lq700(ix,iy) = q700(ixy) ! pwrf_grid%lt700(ix,iy) = t700(ixy) ! pwrf_grid%lu700(ix,iy) = u700(ixy) ! pwrf_grid%lv700(ix,iy) = v700(ixy) ! pwrf_grid%lw700(ix,iy) = w700(ixy) ! pwrf_grid%lz700(ix,iy) = z700(ixy) ! ELSE IF (TRIM(bb2) == "500") THEN ! pwrf_grid%lq500(ix,iy) = q500(ixy) ! pwrf_grid%lt500(ix,iy) = t500(ixy) ! pwrf_grid%lu500(ix,iy) = u500(ixy) ! pwrf_grid%lv500(ix,iy) = v500(ixy) ! pwrf_grid%lw500(ix,iy) = w500(ixy) ! pwrf_grid%lz500(ix,iy) = z500(ixy) ! ELSE IF (TRIM(bb2) == "200")THEN ! pwrf_grid%lq200(ix,iy) = q200(ixy) ! pwrf_grid%lt200(ix,iy) = t200(ixy) ! pwrf_grid%lu200(ix,iy) = u200(ixy) ! pwrf_grid%lv200(ix,iy) = v200(ixy) ! pwrf_grid%lw200(ix,iy) = w200(ixy) ! pwrf_grid%lz200(ix,iy) = z200(ixy) ! ELSE IF (TRIM(bb2) == "100") THEN ! pwrf_grid%lq100(ix,iy) = q100(ixy) ! pwrf_grid%lt100(ix,iy) = t100(ixy) ! pwrf_grid%lu100(ix,iy) = u100(ixy) ! pwrf_grid%lv100(ix,iy) = v100(ixy) ! pwrf_grid%lw100(ix,iy) = w100(ixy) ! pwrf_grid%lz100(ix,iy) = z100(ixy) ! ELSE IF (TRIM(bb2) == "50") THEN ! pwrf_grid%lq50(ix,iy) = q50(ixy) ! pwrf_grid%lt50(ix,iy) = t50(ixy) ! pwrf_grid%lu50(ix,iy) = u50(ixy) ! pwrf_grid%lv50(ix,iy) = v50(ixy) ! pwrf_grid%lw50(ix,iy) = w50(ixy) ! pwrf_grid%lz50(ix,iy) = z50(ixy) ! ELSE IF (TRIM(bb2) == "10") THEN ! pwrf_grid%lq10(ix,iy) = q10(ixy) ! pwrf_grid%lt10(ix,iy) = t10(ixy) ! pwrf_grid%lu10(ix,iy) = u10(ixy) ! pwrf_grid%lv10(ix,iy) = v10(ixy) ! pwrf_grid%lw10(ix,iy) = w10(ixy) ! pwrf_grid%lz10(ix,iy) = z10(ixy) ! END IF ! END DO END DO END DO dx=ddx+2*bdy dy=ddy+2*bdy CALL vect_mat_zstagg(fraca, dx, dy, ddz+1, bdy, pwrf_grid%la_th) CALL vect_mat(beta_prec, dx, dy, ddz, bdy, pwrf_grid%lbeta_prec) var3d=upwd + dnwd+ dnwd0 CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldmc) CALL vect_mat(dnwd, dx, dy, ddz, bdy, pwrf_grid%ldnwd) CALL vect_mat(dnwd0, dx, dy, ddz, bdy, pwrf_grid%ldnwd0) ! Originally is var3d = d_q_ajsb(1:ddx*ddy,1:ddz)/pdtphys var3d = d_q_ajsb/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqajs) var3d = d_q_con/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqcon) CALL vect_mat(d_q_dyn, dx, dy, ddz, bdy, pwrf_grid%ldqdyn) var3d = d_q_eva/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqeva) var3d = d_q_lsc/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlsc) var3d = d_q_lscst/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscst) var3d = d_q_lscth/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqlscth) ! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq! ! MOISTtend in the WRF_LMDZ coupling CALL vect_mat(d_qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%ldqphy) var3d = d_q_ajs/pdtphys - d_q_ajsb/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqthe) var3d = d_q_vdf/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqvdf) var3d = d_q_wake/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldqwak) var3d = d_t_ajsb/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtajs) var3d = d_t_con/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtcon) var3d = d_t_diss/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtdis) CALL vect_mat(d_t_dyn, dx, dy, ddz, bdy, pwrf_grid%ldtdyn) var3d = d_t_ec/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtec) var3d = d_t_eva/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldteva) CALL vect_mat(detr_therm, dx, dy, ddz, bdy, pwrf_grid%ld_th) CALL vect_mat(cldemi, dx, dy, ddz, bdy, pwrf_grid%lcldemi) CALL vect_mat(cldtau, dx, dy, ddz, bdy, pwrf_grid%lcldtau) CALL vect_mat(clwcon, dx, dy, ddz, bdy, pwrf_grid%lclwcon) var3d = d_t_lif/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlif) var3d = d_t_lsc/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlsc) var3d = (d_t_lsc+d_t_eva)/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlschr) var3d = d_t_lscst/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscst) var3d = d_t_lscth/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlscth) var3d=-1.*cool0/RDAY CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlw0) var3d = -1.*cool/RDAY CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtlwr) var3d=d_t_oro/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtoro) CALL vect_mat(d_t, dx, dy, ddz, bdy, pwrf_grid%ldtphy) var3d=heat0/RDAY CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtsw0) var3d=heat/RDAY CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtswr) var3d = d_t_ajs/pdtphys - d_t_ajsb/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtthe) var3d=d_t_vdf/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtvdf) var3d=d_t_wake/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldtwak) var3d=d_u_con/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lducon) CALL vect_mat(d_u_dyn, dx, dy, ddz, bdy, pwrf_grid%ldudyn) var3d=d_u_lif/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldulif) var3d=d_u_oro/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduoro) var3d=d_u_vdf/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lduvdf) var3d=d_v_con/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvcon) CALL vect_mat(d_v_dyn, dx, dy, ddz, bdy, pwrf_grid%ldvdyn) var3d=d_v_lif/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvlif) var3d=d_v_oro/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvoro) var3d=d_v_vdf/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ldvvdf) !NOTfound! CALL vect_mat(ec550aer, dx, dy, ddz, bdy, pwrf_grid%lec550aer) CALL vect_mat(entr_therm, dx, dy, ddz, bdy, pwrf_grid%le_th) CALL vect_mat(coefm(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%levu) CALL vect_mat(fl, dx, dy, ddz, bdy, pwrf_grid%lfl) CALL vect_mat(zphi, dx, dy, ddz, bdy, pwrf_grid%lgeop) var3d = q_seri + ql_seri CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lh2o) CALL vect_mat(fiwc, dx, dy, ddz, bdy, pwrf_grid%liwcon) CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz) CALL vect_mat(coefh(:,:,is_ave), dx, dy, ddz+1, bdy, pwrf_grid%lkz_max) ! CALL vect_mat(lambda_th, dx, dy, ddz, bdy, pwrf_grid%llambda_th) CALL vect_mat(flwc, dx, dy, ddz, bdy, pwrf_grid%llwcon) !NOTknown! CALL vect_mat(ma, dx, dy, ddz, bdy, pwrf_grid%lma) CALL vect_mat(zmasse, dx, dy, ddz, bdy, pwrf_grid%lmass) !NOTknown! CALL vect_mat(mc, dx, dy, ddz, bdy, pwrf_grid%lmc) !NOTknown! CALL vect_mat(mcd, dx, dy, ddz, bdy, pwrf_grid%lmcd) CALL vect_mat(ql_seri, dx, dy, ddz, bdy, pwrf_grid%loliq) CALL vect_mat(q_seri, dx, dy, ddz, bdy, pwrf_grid%lovap) ! L. Fita, LMD January 2014. d_qx is the tendency of moisture due to physiq! ! MOISTtend in the WRF_LMDZ coupling CALL vect_mat(qx(:,:,ivap), dx, dy, ddz, bdy, pwrf_grid%lovapinit) !Notnow! var3d = wo(:,:,1)*dobson_u*1e3 / zmasse / rmo3 * rmd) !Notnow! CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lozone) CALL vect_mat_zstagg(paprs, dx, dy, ddz+1, bdy, pwrf_grid%lpaprs) !NOTknown! CALL vect_mat(pb, dx, dy, ddz, bdy, pwrf_grid%lpb) CALL vect_mat_zstagg(pmflxs, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_i) CALL vect_mat_zstagg(pmflxr, dx, dy, ddz+1, bdy, pwrf_grid%lpr_con_l) CALL vect_mat(pplay, dx, dy, ddz, bdy, pwrf_grid%lpres) CALL vect_mat_zstagg(psfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_i) CALL vect_mat_zstagg(prfl, dx, dy, ddz+1, bdy, pwrf_grid%lpr_lsc_l) var3d = 0. WHERE (ptconv) var3d = 1. CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%lptconv) CALL vect_mat(zqasc, dx, dy, ddz, bdy, pwrf_grid%lq_th) CALL vect_mat(ratqs, dx, dy, ddz, bdy, pwrf_grid%lratqs) CALL vect_mat(re, dx, dy, ddz, bdy, pwrf_grid%lre) CALL vect_mat(ref_ice, dx, dy, ddz, bdy, pwrf_grid%lref_ice) CALL vect_mat(ref_liq, dx, dy, ddz, bdy, pwrf_grid%lref_liq) CALL vect_mat(zx_rh, dx, dy, ddz, bdy, pwrf_grid%lrhum) CALL vect_mat(lwdn, dx, dy, ddz, bdy, pwrf_grid%lrld) CALL vect_mat(lwdn0, dx, dy, ddz, bdy, pwrf_grid%lrldcs) CALL vect_mat(lwup, dx, dy, ddz, bdy, pwrf_grid%lrlu) CALL vect_mat(lwup0, dx, dy, ddz, bdy, pwrf_grid%lrlucs) !NOTknown! CALL vect_mat(rn, dx, dy, ddz, bdy, pwrf_grid%lrn) CALL vect_mat(cldfra, dx, dy, ddz, bdy, pwrf_grid%lrneb) CALL vect_mat(rnebcon, dx, dy, ddz, bdy, pwrf_grid%lrnebcon) CALL vect_mat(rneb, dx, dy, ddz, bdy, pwrf_grid%lrnebls) CALL vect_mat(swdn, dx, dy, ddz, bdy, pwrf_grid%lrsd) CALL vect_mat(swdn0, dx, dy, ddz, bdy, pwrf_grid%lrsdcs) CALL vect_mat(swup, dx, dy, ddz, bdy, pwrf_grid%lrsu) CALL vect_mat(swup0, dx, dy, ddz, bdy, pwrf_grid%lrsucs) CALL vect_mat(fluxt(:,:,1), dx, dy, ddz, bdy, pwrf_grid%lsens_ter) CALL vect_mat(fluxt(:,:,2), dx, dy, ddz, bdy, pwrf_grid%lsens_lic) CALL vect_mat(fluxt(:,:,3), dx, dy, ddz, bdy, pwrf_grid%lsens_oce) CALL vect_mat(fluxt(:,:,4), dx, dy, ddz, bdy, pwrf_grid%lsens_sic) CALL vect_mat(t_seri, dx, dy, ddz, bdy, pwrf_grid%ltemp) CALL vect_mat(theta, dx, dy, ddz, bdy, pwrf_grid%ltheta) var3d=0. DO nsrf=1,nbsrf DO iz=1,ddz+1 var3dstagg(:,iz)=var3dstagg(:,iz)+pctsrf(:,nsrf)*pbl_tke(:,iz,nsrf) ENDDO ENDDO CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke) CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%ltke_max) CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter) CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic) CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce) CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic) !NOTinRegisty! CALL vect_mat(pbl_tke(:,:,1), dx, dy, ddz+1, bdy, pwrf_grid%ltke_ter_max) !NOTinRegisty! CALL vect_mat(pbl_tke(:,:,2), dx, dy, ddz+1, bdy, pwrf_grid%ltke_lic_max) !NOTinRegisty! CALL vect_mat(pbl_tke(:,:,3), dx, dy, ddz+1, bdy, pwrf_grid%ltke_oce_max) !NOTinRegisty! CALL vect_mat(pbl_tke(:,:,4), dx, dy, ddz+1, bdy, pwrf_grid%ltke_sic_max) var3d = d_qx(:,:,ivap)+d_q_dyn CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhus) IF (iflag_thermals == 1) THEN var3d = d_q_con/pdtphys ELSE IF ((iflag_thermals > 1) .AND. (iflag_wake == 1)) THEN var3d = d_q_con/pdtphys + d_q_ajs/pdtphys + d_q_wake/pdtphys END IF CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusc) var3d = d_q_lsc/pdtphys + d_q_eva/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnhusscpbl) var3d = d_t + d_t_dyn CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltnt) IF (iflag_thermals == 1) THEN var3d = d_t_con/pdtphys + d_t_ajsb/pdtphys ELSE IF ((iflag_thermals > 1) .AND. (iflag_wake == 1)) THEN var3d=d_t_con/pdtphys + d_t_ajs/pdtphys + d_t_wake/pdtphys END IF CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntc) var3d = heat/RDAY - cool/RDAY CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntr) var3d = (d_t_lsc + d_t_eva + d_t_vdf)/pdtphys CALL vect_mat(var3d, dx, dy, ddz, bdy, pwrf_grid%ltntscpbl) CALL vect_mat(upwd, dx, dy, ddz, bdy, pwrf_grid%lupwd) CALL vect_mat(u_seri, dx, dy, ddz, bdy, pwrf_grid%lvitu) CALL vect_mat(v_seri, dx, dy, ddz, bdy, pwrf_grid%lvitv) CALL vect_mat(omega, dx, dy, ddz, bdy, pwrf_grid%lvitw) CALL vect_mat_zstagg(Vprecip, dx, dy, ddz+1, bdy, pwrf_grid%lvprecip) CALL vect_mat(wake_deltaq, dx, dy, ddz, bdy, pwrf_grid%lwake_deltaq) CALL vect_mat(wake_deltat, dx, dy, ddz, bdy, pwrf_grid%lwake_deltat) CALL vect_mat(wake_omg, dx, dy, ddz, bdy, pwrf_grid%lwake_omg) CALL vect_mat(wdtrainA, dx, dy, ddz, bdy, pwrf_grid%lwdtraina) CALL vect_mat(wdtrainM, dx, dy, ddz, bdy, pwrf_grid%lwdtrainm) ! L. Fita, LMD, January 2014. pphis is the given value of the pressure at the surface var3dstagg(:,1) = pphis(:)/RG DO iz=1,ddz var3dstagg(:,iz+1)= var3dstagg(:,iz)-(t_seri(:,iz)*RD*(paprs(:,iz+1)- & paprs(:,iz)) )/ (pplay(:,iz)*RG) ENDDO CALL vect_mat(var3dstagg, dx, dy, ddz, bdy, pwrf_grid%lzfull) var3dstagg(:,1) = pphis(:)/RG-( (t_seri(:,1)+zxtsol)/2.*RD*(pplay(:,1)- & paprs(:,1)) )/ ((paprs(:,1)+pplay(:,1))/2.*RG) DO iz=1, ddz-1 var3dstagg(:,iz+1)= var3dstagg(:,iz)-( (t_seri(:,iz)+t_seri(:,iz+1))/ & 2.*RD*(pplay(:,iz+1)-pplay(:,iz)))/( paprs(:,iz)*RG) ENDDO CALL vect_mat_zstagg(var3dstagg, dx, dy, ddz+1, bdy, pwrf_grid%lzhalf) RETURN END SUBROUTINE get_lmdz_out SUBROUTINE put_lmdz_in_WRFout(is, ie, js, je, ks, ke, ddx, ddy, ddz, bdy, dlmdz,xp,& yp, lmdp, Nks, Nz0, nqtot, ivap, iliq, iflag_pbl, iflag_con, iflag_thermals, & iflag_wake, itap, pdtphys, paprs, pplay, pphi, pphis, presnivs, qx, d_u, d_v, & d_t, d_qx, d_ps, orchidlev, tsoilorchid, nslayers, wslay, wzs, wdzs, wnest_pos, & wq2, wt2, wth2, wpsfc, wu10, wv10, wresm, wzetatop, wtslb, wsmois, & wsh2o, wsmcrel, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh, & wcanwat, wsst, wsstsk, wlai, wh_diabatic, wsina, wtsk, wtiso, & wmax_msftx, wmax_msfty, wrainc, wraincv, wrainsh, wrainshv, wrainnc, wrainncv, & wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc, whailncv, wstepave_count, & wcldfra, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx, & wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb) ! wsave_topo_from_real, wavg_fuel_frac, wuah, wvah, wseed1, wseed2) ! Subroutine to get LMDZ output and get into WRF standard output variables ! (as they appear in the wrfout) ! LMDZ modules USE indice_sol_mod USE phys_state_var_mod USE phys_output_var_mod ! USE pbl_surface_mod USE phys_local_var_mod USE comgeomphy ! WRF_LMDZ modules USE output_lmdz_NOmodule USE lmdz_wrf_variables_mod IMPLICIT NONE INTEGER, INTENT(IN) :: is, ie, js, je, ks, ke, & ddx, ddy, ddz, bdy, dlmdz, xp, yp, lmdp, Nks, Nz0, nqtot, ivap, iliq, & iflag_pbl, iflag_con, iflag_thermals, iflag_wake, itap, nslayers REAL, INTENT(IN) :: pdtphys REAL, DIMENSION(dlmdz), INTENT(IN) :: pphis, d_ps REAL, DIMENSION(ddz), INTENT(IN) :: presnivs REAL, DIMENSION(dlmdz,ddz+1), INTENT(IN) :: paprs REAL, DIMENSION(dlmdz,ddz), INTENT(IN) :: pplay, pphi, d_u, d_v, d_t REAL, DIMENSION(dlmdz,ddz,nqtot), INTENT(IN) :: qx, d_qx REAL, DIMENSION(Nz0), INTENT(IN) :: orchidlev REAL, DIMENSION(dlmdz,Nz0,Nks), INTENT(IN) :: tsoilorchid INTEGER, INTENT(INOUT) :: wstepave_count REAL, INTENT(INOUT) :: wresm, wzetatop, wtiso, & wmax_msftx, wmax_msfty REAL, DIMENSION(1:nslayers), INTENT(INOUT) :: wzs, wdzs REAL, DIMENSION(1:nslayers), INTENT(IN) :: wslay REAL, DIMENSION(is:ie,js:je), INTENT(INOUT) :: wnest_pos, wq2, wt2, & wth2, wpsfc, wu10, wv10, wsfroff, wudroff, wgrdflx, wacgrdflx, wsnow, wsnowh, & wcanwat, wsst, wsstsk, wlai, wsina, wtsk, wrainc, wraincv, wrainsh, wrainshv, & wrainnc, wrainncv, wsnownc, wsnowncv, wgraupelnc, wgraupelncv, whailnc, & whailncv, wswdown, wglw, wswnorm, wolr, wemiss, wtmn, wust, wpblh, whfx, & wqfx, wlh, wachfx, waclhf, wsnowc, wsr, wpotevp, wsnopcx, wsoiltb REAL, DIMENSION(is:ie,1:nslayers,js:je), INTENT(INOUT) :: wtslb, wsmois, wsh2o, & wsmcrel REAL, DIMENSION(is:ie,ks:ke,js:je), INTENT(INOUT) :: wh_diabatic, wcldfra ! EXTERNAL suphel ! Local INTEGER :: dx, dy, ix, iy, ixy, iz INTEGER :: nsrf REAL :: varixy REAL, DIMENSION(Nz0) :: yvals1 REAL :: RDAY, RG, RD CHARACTER(LEN=50) :: fname !!!!!!! Variables ! pwrf_grid: WRF grid structure fname = "'put_lmdz_in_WRFout'" ! L. Fita, LMD January 2014. It should work using suphel, but not? ! CALL suphel RDAY = 86400. RG = 9.81 RD = 287. wdzs(1)=wslay(nslayers) DO iz=1,nslayers-1 wdzs(iz+1) = wslay(nslayers-iz+1) - wslay(nslayers-iz) END DO ! wmax_msftx = ! wmax_msfty = ! wresm = ! wstepave_count = ! wtiso = ! wzetatop = wzs = wslay DO ix=1,ddx DO iy=1,ddy ixy=ddx*(iy-1)+ix wacgrdflx(ix,iy) = wacgrdflx(ix,iy) + 0. wachfx(ix,iy) = wachfx(ix,iy) + sens(ixy) waclhf(ix,iy) = waclhf(ix,iy) + zxfluxlat(ixy) DO iz=1,ddz wcldfra(ix,iz,iy) = cldfra(ixy,iz) END DO wcanwat(ix,iy) = 0. wglw(ix,iy) = sollwdown(ixy) ! These are precipitation fluxes !! ! wgraupelnc(ix,iy) = wgraupelnc(ix,iy) + prmflxs(ixy) + prfs(ixy,iz) ! wgraupelncv(ix,iy) = prmflxs(ixy) + prfs(ixy,iz) wgrdflx(ix,iy) = 0. whailnc(ix,iy) = whailnc(ix,iy) + 0. whailncv(ix,iy) = 0. wh_diabatic(ix,ks:ke,iy) = 0. whfx(ix,iy) = sens(ixy) wlai(ix,iy) = 0. wlh(ix,iy) = zxfluxlat(ixy) wnest_pos(ix,iy) = 0. !! wnoahres(ix,iy) = 0. wolr(ix,iy) = toplw(ixy) wpblh(ix,iy) = s_pblh(ixy) ! Should be something from evappot_srf ? wpotevp(ix,iy) = wpotevp(ix,iy) + 0. wpsfc(ix,iy) = paprs(ixy,1) wq2(ix,iy) = zq2m(ixy) wqfx(ix,iy) = 0. wrainc(ix,iy) = wrainc(ix,iy) + (rain_con(ixy) + snow_con(ixy))*pdtphys wraincv(ix,iy) = (rain_con(ixy) + snow_con(ixy))*pdtphys wrainnc(ix,iy) = wrainnc(ix,iy) + (rain_lsc(ixy)+snow_lsc(ixy))*pdtphys wrainncv(ix,iy) = (rain_lsc(ixy)+snow_lsc(ixy))*pdtphys wrainsh(ix,iy) = wrainsh(ix,iy) + 0. wrainshv(ix,iy) = 0. ! L. Fita, LMD. February 2014 ! LMDZ run_off_lic is for continental ice. We could try to get run_off_ter, but.... ! wsfroff(ix,iy) = 0. DO iz=1,Nks yvals1(iz) = tsoilorchid(ixy,iz,1)*pctsrf(ixy,is_ter) + & tsoilorchid(ixy,iz,2)*pctsrf(ixy,is_lic) + & tsoilorchid(ixy,iz,3)*pctsrf(ixy,is_oce) + & tsoilorchid(ixy,iz,4)*pctsrf(ixy,is_sic) END DO CALL interpolate_layers(Nz0,orchidlev,yvals1,nslayers,wslay,wtslb(ix,:,iy)) ! No way ??? ! wsh2o(is:ie,iz,js:je) = wsh2o(is:ie,iz,js:je) + ! wsmcrel(is:ie,iz,js:je) = 0. ! wsinalpha(ix,iy) = wsina(ix,iy) = 0. wsnopcx(ix,iy) = wsnopcx(ix,iy) + 0. !! wsnowh(ix,iy) = snowght(ixy) wsnownc(ix,iy) = wsnownc(ix,iy) + snow_fall(ixy)*pdtphys wsnowncv(ix,iy) = snow_fall(ixy)*pdtphys wsoiltb(ix,iy) = wtslb(ix,nslayers,iy) IF (rain_fall(ixy) /= 0.) THEN wsr(ix,iy) = ( snow_fall(ixy) ) / ( rain_fall(ixy)) ELSE wsr(ix,iy) = 0. END IF wsstsk(ix,iy) = ftsol(ixy,is_oce) wsst(ix,iy) = sst(ixy) wswdown(ix,iy) = swdn(ixy,1)*pctsrf(ixy,is_ter) + & swdn(ixy,2)*pctsrf(ixy,is_lic) + swdn(ixy,3)*pctsrf(ixy,is_oce) + & swdn(ixy,4)*pctsrf(ixy,is_sic) wswnorm(ix,iy) = 0. wt2(ix,iy) = zt2m(ixy) ! This is not exact!! CALL temp_to_theta1(zt2m(ixy), paprs(ixy,1), wth2(ix,iy)) wtsk(ix,iy) = ftsol(ixy, is_ter)*pctsrf(ixy,is_ter) + ftsol(ixy, is_lic)* & pctsrf(ixy,is_lic) + ftsol(ixy, is_oce)*pctsrf(ixy,is_oce) + & ftsol(ixy, is_sic)*pctsrf(ixy,is_sic) wu10(ix,iy) = zu10m(ixy) wudroff(ix,iy) = wudroff(ix,iy) + 0. wv10(ix,iy) = zv10m(ixy) END DO END DO END SUBROUTINE put_lmdz_in_WRFout END MODULE wrf_lmdz_mod