Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/etat0phys_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/dynlonlat_phylonlat/etat0phys_netcdf.F90 (revision 2300) +++ LMDZ5/trunk/libf/dynlonlat_phylonlat/etat0phys_netcdf.F90 (revision 2300) @@ -0,0 +1,558 @@ +MODULE etat0phys +! +!******************************************************************************* +! Purpose: Create physical initial state using atmospheric fields from a +! database of atmospheric to initialize the model. +!------------------------------------------------------------------------------- +! Comments: +! +! * This module is designed to work for Earth (and with ioipsl) +! +! * etat0phys_netcdf routine can access to NetCDF data through subroutines: +! "start_init_phys" for variables contained in file "ECPHY.nc": +! 'ST' : Surface temperature +! 'CDSW' : Soil moisture +! "start_init_orog" for variables contained in file "Relief.nc": +! 'RELIEF' : High resolution orography +! +! * The land mask and corresponding weights can be: +! 1) computed using the ocean mask from the ocean model (to ensure ocean +! fractions are the same for atmosphere and ocean) for coupled runs. +! File name: "o2a.nc" ; variable name: "OceMask" +! 2) computed from topography file "Relief.nc" for forced runs. +! +! * Allowed values for read_climoz flag are 0, 1 and 2: +! 0: do not read an ozone climatology +! 1: read a single ozone climatology that will be used day and night +! 2: read two ozone climatologies, the average day and night climatology +! and the daylight climatology +!------------------------------------------------------------------------------- +! * There is a big mess with the longitude size. Should it be iml or iml+1 ? +! I have chosen to use the iml+1 as an argument to this routine and we declare +! internaly smaller fields when needed. This needs to be cleared once and for +! all in LMDZ. A convention is required. +!------------------------------------------------------------------------------- + + USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo + USE assert_eq_m, ONLY: assert_eq +#ifdef CPP_EARTH + USE dimphy + USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, & + rlon, solsw, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, & + rlat, sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, & + sig1, ftsol, clwcon, fm_therm, wake_Cstar, pctsrf, entr_therm,radpas, f0,& + zmax0,fevap, rnebcon,falb_dir, wake_fip, agesno, detr_therm, pbl_tke, & + phys_state_var_init +#endif + + PRIVATE + PUBLIC :: etat0phys_netcdf + + include "iniprint.h" + include "dimensions.h" + include "paramet.h" + include "comgeom2.h" + include "comvert.h" + include "comconst.h" + include "dimsoil.h" + include "temps.h" + include "comdissnew.h" + include "serre.h" + include "clesphys.h" + REAL, SAVE :: deg2rad + REAL, SAVE, ALLOCATABLE :: tsol(:) +! REAL, SAVE, ALLOCATABLE :: rugo(:,:) ! ??? COMPUTED BUT NOT USED ??? + INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys + REAL, ALLOCATABLE, SAVE :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:) + CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc" + CHARACTER(LEN=256), PARAMETER :: title="RELIEF" + + +CONTAINS + + +!------------------------------------------------------------------------------- +! +SUBROUTINE etat0phys_netcdf(ib, masque, phis) +! +!------------------------------------------------------------------------------- +! Purpose: Creates initial states +!------------------------------------------------------------------------------- +! Note: This routine is designed to work for Earth +!------------------------------------------------------------------------------- + USE control_mod +#ifdef CPP_EARTH + USE infotrac + USE fonte_neige_mod + USE pbl_surface_mod + USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz + USE indice_sol_mod + USE conf_phys_m, ONLY: conf_phys + USE exner_hyb_m, ONLY: exner_hyb + USE exner_milieu_m, ONLY: exner_milieu + USE test_disvert_m, ONLY: test_disvert + USE grid_atob_m, ONLY: grille_m +#endif + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation + REAL, INTENT(INOUT) :: masque(:,:) !--- Land mask dim(iip1,jjp1) + REAL, INTENT(INOUT) :: phis (:,:) !--- Ground geopotential dim(iip1,jjp1) +#ifndef CPP_EARTH + WRITE(lunout,*)'etat0phys_netcdf: Earth-specific routine, needs Earth physics' +#else +!------------------------------------------------------------------------------- +! Local variables: + CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt + INTEGER :: i, j, l, ji, iml, jml + REAL :: phystep + REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp + REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder + REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf + REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil + +!--- Local variables for sea-ice reading: + LOGICAL :: read_mask + INTEGER :: iml_lic, jml_lic, isst(klon-2) + INTEGER :: fid, llm_tmp, ttm_tmp, itaul(1) + REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:) + REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:) + REAL :: date, lev(1), dummy + REAL :: flic_tmp(SIZE(masque,1),SIZE(masque,2)) + +!--- Arguments for conf_phys + LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats + REAL :: solarlong0, seuil_inversion, fact_cldcon, facttemps + LOGICAL :: ok_newmicro + INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs + REAL :: ratqsbas, ratqshaut, tau_ratqs + LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple + INTEGER :: flag_aerosol + LOGICAL :: flag_aerosol_strat + LOGICAL :: new_aod + REAL :: bl95_b0, bl95_b1 + INTEGER :: read_climoz !--- Read ozone climatology + REAL :: alp_offset + + deg2rad= pi/180.0 + iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml") + jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml") + +! Grid construction and miscellanous initializations. +!******************************************************************************* + CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & + callstats, & + solarlong0,seuil_inversion, & + fact_cldcon, facttemps,ok_newmicro,iflag_radia, & + iflag_cldcon, & + iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & + ok_ade, ok_aie, ok_cdnc, aerosol_couple, & + flag_aerosol, flag_aerosol_strat, new_aod, & + bl95_b0, bl95_b1, & + read_climoz, & + alp_offset) + + CALL phys_state_var_init(read_climoz) + + co2_ppm0 = co2_ppm !--- Initial atmospheric CO2 conc. from .def file + + rlat(1) = pi/2. + DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO + rlat(klon) = - pi/2. + rlat(:)=rlat(:)*(180.0/pi) + + rlon(1) = 0.0 + DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO + rlon(klon) = 0.0 + rlon(:)=rlon(:)*(180.0/pi) + +! Compute ground geopotential, sub-cells quantities and possibly the mask. +!******************************************************************************* + read_mask=ANY(masque/=-99999.); masque_tmp=masque + CALL start_init_orog(rlonv, rlatu, phis, masque_tmp) + WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt) + IF(.NOT.read_mask) THEN !--- Keep mask form orography + masque=masque_tmp + IF(prt_level>=1) THEN + WRITE(lunout,*)'BUILT MASK :' + WRITE(lunout,fmt) NINT(masque) + END IF + WHERE( masque(:,:)=1) & !--- Ozone climatology + CALL regr_lat_time_climoz(read_climoz) + +! Sub-surfaces initialization +!******************************************************************************* +!--- Read and interpolate on model T-grid soil fraction and soil ice fraction. + CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) + ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) + ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) + ALLOCATE( fraclic(iml_lic,jml_lic)) + CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & + & lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) + CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) + CALL flinclo(fid) + WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic + IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees + IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. + dlon_lic(:)=lon_lic(:,1) + dlat_lic(:)=lat_lic(1,:) + CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) ) + flic_tmp(iml,:)=flic_tmp(1,:) + +!--- To the physical grid + pctsrf(:,:) = 0. + CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic)) + +!--- Adequation with soil/sea mask + WHERE(pctsrf(:,is_lic)EPSFRA) THEN + IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN + pctsrf(ji,is_lic)=zmasq(ji) + pctsrf(ji,is_ter)=0. + ELSE + pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) + IF(pctsrf(ji,is_ter)0.) isst=1 + +!--- It is checked that the sub-surfaces sum is equal to 1. + ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) + IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points' + +! Write physical initial state +!******************************************************************************* + WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad + phystep = dtvr * FLOAT(iphysiq) + radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) + WRITE(lunout,*)'phystep =', phystep, radpas + +! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0 +!******************************************************************************* + DO i=1,nbsrf; ftsol(:,i) = tsol; END DO + DO i=1,nbsrf; snsrf(:,i) = sn; END DO + falb_dir(:,is_ter,:) = 0.08 + falb_dir(:,is_lic,:) = 0.6 + falb_dir(:,is_oce,:) = 0.5 + falb_dir(:,is_sic,:) = 0.6 + fevap(:,:) = 0. + DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO + DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO + rain_fall = 0. + snow_fall = 0. + solsw = 165. + sollw = -53. + t_ancien = 273.15 + q_ancien = 0. + agesno = 0. + + z0m(:,is_oce) = rugmer(:) + z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) + z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) + z0m(:,is_sic) = 0.001 + z0h(:,:)=z0m(:,:) + + fder = 0.0 + clwcon = 0.0 + rnebcon = 0.0 + ratqs = 0.0 + run_off_lic_0 = 0.0 + rugoro = 0.0 + +! Before phyredem calling, surface modules and values to be saved in startphy.nc +! are initialized +!******************************************************************************* + dummy = 1.0 + pbl_tke(:,:,:) = 1.e-8 + zmax0(:) = 40. + f0(:) = 1.e-5 + sig1(:,:) = 0. + w01(:,:) = 0. + wake_deltat(:,:) = 0. + wake_deltaq(:,:) = 0. + wake_s(:) = 0. + wake_cstar(:) = 0. + wake_fip(:) = 0. + wake_pe = 0. + fm_therm = 0. + entr_therm = 0. + detr_therm = 0. + + CALL fonte_neige_init(run_off_lic_0) + CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil ) + CALL phyredem( "startphy.nc" ) + +! WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' +! WRITE(lunout,*)'entree histclo' + CALL histclo() + +#endif +!#endif of #ifdef CPP_EARTH + +END SUBROUTINE etat0phys_netcdf +! +!------------------------------------------------------------------------------- + + +#ifdef CPP_EARTH +!------------------------------------------------------------------------------- +! +SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque) +! +!=============================================================================== +! Comment: +! This routine launch grid_noro, which computes parameters for SSO scheme as +! described in LOTT & MILLER (1997) and LOTT(1999). +!=============================================================================== + USE conf_dat_m, ONLY: conf_dat2d +! USE grid_atob_m, ONLY: rugsoro + USE grid_noro_m, ONLY: grid_noro + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) + REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) +!------------------------------------------------------------------------------- +! Local variables: + CHARACTER(LEN=256) :: modname="start_init_orog" + CHARACTER(LEN=256) :: title="RELIEF" + INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) + REAL :: lev(1), date, dt + REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) + REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var (:,:) + REAL, ALLOCATABLE :: zmea0(:,:), zstd0(:,:), zsig0(:,:) + REAL, ALLOCATABLE :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:) +!------------------------------------------------------------------------------- + iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") + jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") + +!--- HIGH RESOLUTION OROGRAPHY + CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) + + ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) + CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& + lev, ttm_tmp, itau, date, dt, fid) + ALLOCATE(relief_hi(iml_rel,jml_rel)) + CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) + CALL flinclo(fid) + +!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS + ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) + lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad + lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad + +!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS + ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) + CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) + DEALLOCATE(lon_ini,lat_ini) + +!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro + WRITE(lunout,*) + WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***' + +!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES + ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation + ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy + ALLOCATE(zthe0(iml,jml)) !--- Highest slope orientation + ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights + +!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS + CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0, & + zsig0,zgam0,zthe0,zpic0,zval0,masque) + phis = phis * 9.81 + phis(iml,:) = phis(1,:) + +!--- COMPUTE SURFACE ROUGHNESS +! WRITE(lunout,*) +! WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***' +! ALLOCATE(tmp_var(iml-1,jml)) +! CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var) +! ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:) +! DEALLOCATE(tmp_var) + DEALLOCATE(relief_hi,lon_rad,lat_rad) + +!--- PUT QUANTITIES TO PHYSICAL GRID + CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0) + CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0) + CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0) + CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0) + CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0) + CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0) + CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0) + + +END SUBROUTINE start_init_orog +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis) +! +!=============================================================================== +! Purpose: Compute tsol and qsol, knowing phis. +!=============================================================================== + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) + REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) + LOGICAL, INTENT(IN) :: ibar + REAL, INTENT(IN) :: phis(:,:) ! dim (iml,jml) +!------------------------------------------------------------------------------- +! Local variables: + CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc" + REAL :: date, dt + INTEGER :: iml, jml, jml2, itau(1) + REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) + REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) + REAL, ALLOCATABLE :: ts(:,:), qs(:,:) +!------------------------------------------------------------------------------- + iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml") + jml=assert_eq(SIZE(lat_in),SIZE(phis,2), TRIM(modname)//" jml") + jml2=SIZE(lat_in2) + + WRITE(lunout,*)'Opening the surface analysis' + CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys) + WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys + + ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys)) + ALLOCATE(levphys_ini(llm_phys)) + CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys, & + lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys) + +!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS + ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys)) + lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad + lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad + + ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) + CALL get_var_phys('ST' ,ts) !--- SURFACE TEMPERATURE + CALL get_var_phys('CDSW',qs) !--- SOIL MOISTURE + CALL flinclo(fid_phys) + DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) + +!--- TSOL AND QSOL ON PHYSICAL GRID + ALLOCATE(tsol(klon)) + CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol) + CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol) + DEALLOCATE(ts,qs) + +CONTAINS + +!------------------------------------------------------------------------------- +! +SUBROUTINE get_var_phys(title,field) +! +!------------------------------------------------------------------------------- + USE conf_dat_m, ONLY: conf_dat2d + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + CHARACTER(LEN=*), INTENT(IN) :: title + REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: tllm +!------------------------------------------------------------------------------- + SELECT CASE(title) + CASE('SP'); tllm=0 + CASE('ST','CDSW'); tllm=llm_phys + END SELECT + IF(ALLOCATED(field)) RETURN + ALLOCATE(field(iml,jml)); field(:,:)=0. + CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana) + CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar) + CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana, & + lon_in, lat_in, lon_in2, lat_in2, field) + +END SUBROUTINE get_var_phys +! +!------------------------------------------------------------------------------- +! +END SUBROUTINE start_init_phys +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) +! +!------------------------------------------------------------------------------- + USE inter_barxy_m, ONLY: inter_barxy + USE grid_atob_m, ONLY: grille_m + IMPLICIT NONE +!------------------------------------------------------------------------------- +! Arguments: + CHARACTER(LEN=*), INTENT(IN) :: nam + LOGICAL, INTENT(IN) :: ibar, ibeg + REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) + REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) + REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) + REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) + REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) +!------------------------------------------------------------------------------- +! Local variables: + CHARACTER(LEN=256) :: modname="interp_startvar" + INTEGER :: ii, jj, i1, j1, j2 + REAL, ALLOCATABLE :: vtmp(:,:) +!------------------------------------------------------------------------------- + ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") + jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") + i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") + j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") + j2=SIZE(lat2) + ALLOCATE(vtmp(i1-1,j1)) + IF(ibar) THEN + IF(ibeg.AND.prt_level>1) THEN + WRITE(lunout,*)"--------------------------------------------------------" + WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" + WRITE(lunout,*)"--------------------------------------------------------" + END IF + CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) + ELSE + CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) + END IF + CALL gr_int_dyn(vtmp, varo, i1-1, j1) + +END SUBROUTINE interp_startvar +! +!------------------------------------------------------------------------------- + +#endif +!#endif of #ifdef CPP_EARTH + +END MODULE etat0phys +! +!******************************************************************************* + Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/ioipsl_getin_p_mod.F90 =================================================================== --- LMDZ5/trunk/libf/dynlonlat_phylonlat/ioipsl_getin_p_mod.F90 (revision 2299) +++ (revision ) @@ -1,169 +1,0 @@ -! -! $Id$ -! -MODULE ioipsl_getin_p_mod -! To use getin in a parallel context -!--------------------------------------------------------------------- -#ifdef CPP_IOIPSL -USE ioipsl, ONLY: getin -#else -USE ioipsl_getincom, ONLY: getin -#endif -USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root -USE mod_phys_lmdz_omp_data, ONLY : is_omp_root -use mod_phys_lmdz_para, only : bcast -!- -IMPLICIT NONE -!- -PRIVATE -PUBLIC :: getin_p -!- -INTERFACE getin_p - - MODULE PROCEDURE getinrs_p, getinr1d_p, getinr2d_p, & - & getinis_p, getini1d_p, getini2d_p, & - & getincs_p, & - & getinls_p, getinl1d_p, getinl2d_p -END INTERFACE -!- -CONTAINS - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!! Definition des getin -> bcast !! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!! -- Les chaines de caracteres -- !! - - SUBROUTINE getincs_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - CHARACTER(LEN=*),INTENT(INOUT) :: VarOut - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getincs_p - -!! -- Les entiers -- !! - - SUBROUTINE getinis_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - INTEGER,INTENT(INOUT) :: VarOut - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinis_p - - SUBROUTINE getini1d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - INTEGER,INTENT(INOUT) :: VarOut(:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getini1d_p - - SUBROUTINE getini2d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - INTEGER,INTENT(INOUT) :: VarOut(:,:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getini2d_p - -!! -- Les flottants -- !! - - SUBROUTINE getinrs_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - REAL,INTENT(INOUT) :: VarOut - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinrs_p - - SUBROUTINE getinr1d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - REAL,INTENT(INOUT) :: VarOut(:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinr1d_p - - SUBROUTINE getinr2d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - REAL,INTENT(INOUT) :: VarOut(:,:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinr2d_p - -!! -- Les Booleens -- !! - - SUBROUTINE getinls_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - LOGICAL,INTENT(INOUT) :: VarOut - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinls_p - - SUBROUTINE getinl1d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - LOGICAL,INTENT(INOUT) :: VarOut(:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinl1d_p - - SUBROUTINE getinl2d_p(VarIn,VarOut) - IMPLICIT NONE - CHARACTER(LEN=*),INTENT(IN) :: VarIn - LOGICAL,INTENT(INOUT) :: VarOut(:,:) - -!$OMP BARRIER - IF (is_mpi_root .AND. is_omp_root) THEN - CALL getin(VarIn,VarOut) - ENDIF - CALL bcast(VarOut) - END SUBROUTINE getinl2d_p -!- -!----------------------------- -!----------------------------- -!----------------------------- - -END MODULE ioipsl_getin_p_mod - Index: LMDZ5/trunk/libf/misc/ioipsl_getin_p_mod.F90 =================================================================== --- LMDZ5/trunk/libf/misc/ioipsl_getin_p_mod.F90 (revision 2300) +++ LMDZ5/trunk/libf/misc/ioipsl_getin_p_mod.F90 (revision 2300) @@ -0,0 +1,169 @@ +! +! $Id$ +! +MODULE ioipsl_getin_p_mod +! To use getin in a parallel context +!--------------------------------------------------------------------- +#ifdef CPP_IOIPSL +USE ioipsl, ONLY: getin +#else +USE ioipsl_getincom, ONLY: getin +#endif +USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root +USE mod_phys_lmdz_omp_data, ONLY : is_omp_root +use mod_phys_lmdz_para, only : bcast +!- +IMPLICIT NONE +!- +PRIVATE +PUBLIC :: getin_p +!- +INTERFACE getin_p + + MODULE PROCEDURE getinrs_p, getinr1d_p, getinr2d_p, & + & getinis_p, getini1d_p, getini2d_p, & + & getincs_p, & + & getinls_p, getinl1d_p, getinl2d_p +END INTERFACE +!- +CONTAINS + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Definition des getin -> bcast !! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!! -- Les chaines de caracteres -- !! + + SUBROUTINE getincs_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + CHARACTER(LEN=*),INTENT(INOUT) :: VarOut + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getincs_p + +!! -- Les entiers -- !! + + SUBROUTINE getinis_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + INTEGER,INTENT(INOUT) :: VarOut + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinis_p + + SUBROUTINE getini1d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + INTEGER,INTENT(INOUT) :: VarOut(:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getini1d_p + + SUBROUTINE getini2d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + INTEGER,INTENT(INOUT) :: VarOut(:,:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getini2d_p + +!! -- Les flottants -- !! + + SUBROUTINE getinrs_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + REAL,INTENT(INOUT) :: VarOut + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinrs_p + + SUBROUTINE getinr1d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + REAL,INTENT(INOUT) :: VarOut(:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinr1d_p + + SUBROUTINE getinr2d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + REAL,INTENT(INOUT) :: VarOut(:,:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinr2d_p + +!! -- Les Booleens -- !! + + SUBROUTINE getinls_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + LOGICAL,INTENT(INOUT) :: VarOut + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinls_p + + SUBROUTINE getinl1d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + LOGICAL,INTENT(INOUT) :: VarOut(:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinl1d_p + + SUBROUTINE getinl2d_p(VarIn,VarOut) + IMPLICIT NONE + CHARACTER(LEN=*),INTENT(IN) :: VarIn + LOGICAL,INTENT(INOUT) :: VarOut(:,:) + +!$OMP BARRIER + IF (is_mpi_root .AND. is_omp_root) THEN + CALL getin(VarIn,VarOut) + ENDIF + CALL bcast(VarOut) + END SUBROUTINE getinl2d_p +!- +!----------------------------- +!----------------------------- +!----------------------------- + +END MODULE ioipsl_getin_p_mod + Index: LMDZ5/trunk/libf/phylmd/etat0phys_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/etat0phys_netcdf.F90 (revision 2299) +++ (revision ) @@ -1,558 +1,0 @@ -MODULE etat0phys -! -!******************************************************************************* -! Purpose: Create physical initial state using atmospheric fields from a -! database of atmospheric to initialize the model. -!------------------------------------------------------------------------------- -! Comments: -! -! * This module is designed to work for Earth (and with ioipsl) -! -! * etat0phys_netcdf routine can access to NetCDF data through subroutines: -! "start_init_phys" for variables contained in file "ECPHY.nc": -! 'ST' : Surface temperature -! 'CDSW' : Soil moisture -! "start_init_orog" for variables contained in file "Relief.nc": -! 'RELIEF' : High resolution orography -! -! * The land mask and corresponding weights can be: -! 1) computed using the ocean mask from the ocean model (to ensure ocean -! fractions are the same for atmosphere and ocean) for coupled runs. -! File name: "o2a.nc" ; variable name: "OceMask" -! 2) computed from topography file "Relief.nc" for forced runs. -! -! * Allowed values for read_climoz flag are 0, 1 and 2: -! 0: do not read an ozone climatology -! 1: read a single ozone climatology that will be used day and night -! 2: read two ozone climatologies, the average day and night climatology -! and the daylight climatology -!------------------------------------------------------------------------------- -! * There is a big mess with the longitude size. Should it be iml or iml+1 ? -! I have chosen to use the iml+1 as an argument to this routine and we declare -! internaly smaller fields when needed. This needs to be cleared once and for -! all in LMDZ. A convention is required. -!------------------------------------------------------------------------------- - - USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo - USE assert_eq_m, ONLY: assert_eq -#ifdef CPP_EARTH - USE dimphy - USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, & - rlon, solsw, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, & - rlat, sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, & - sig1, ftsol, clwcon, fm_therm, wake_Cstar, pctsrf, entr_therm,radpas, f0,& - zmax0,fevap, rnebcon,falb_dir, wake_fip, agesno, detr_therm, pbl_tke, & - phys_state_var_init -#endif - - PRIVATE - PUBLIC :: etat0phys_netcdf - - include "iniprint.h" - include "dimensions.h" - include "paramet.h" - include "comgeom2.h" - include "comvert.h" - include "comconst.h" - include "dimsoil.h" - include "temps.h" - include "comdissnew.h" - include "serre.h" - include "clesphys.h" - REAL, SAVE :: deg2rad - REAL, SAVE, ALLOCATABLE :: tsol(:) -! REAL, SAVE, ALLOCATABLE :: rugo(:,:) ! ??? COMPUTED BUT NOT USED ??? - INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys - REAL, ALLOCATABLE, SAVE :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:) - CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc" - CHARACTER(LEN=256), PARAMETER :: title="RELIEF" - - -CONTAINS - - -!------------------------------------------------------------------------------- -! -SUBROUTINE etat0phys_netcdf(ib, masque, phis) -! -!------------------------------------------------------------------------------- -! Purpose: Creates initial states -!------------------------------------------------------------------------------- -! Note: This routine is designed to work for Earth -!------------------------------------------------------------------------------- - USE control_mod -#ifdef CPP_EARTH - USE infotrac - USE fonte_neige_mod - USE pbl_surface_mod - USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz - USE indice_sol_mod - USE conf_phys_m, ONLY: conf_phys - USE exner_hyb_m, ONLY: exner_hyb - USE exner_milieu_m, ONLY: exner_milieu - USE test_disvert_m, ONLY: test_disvert - USE grid_atob_m, ONLY: grille_m -#endif - IMPLICIT NONE -!------------------------------------------------------------------------------- -! Arguments: - LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation - REAL, INTENT(INOUT) :: masque(:,:) !--- Land mask dim(iip1,jjp1) - REAL, INTENT(INOUT) :: phis (:,:) !--- Ground geopotential dim(iip1,jjp1) -#ifndef CPP_EARTH - WRITE(lunout,*)'etat0phys_netcdf: Earth-specific routine, needs Earth physics' -#else -!------------------------------------------------------------------------------- -! Local variables: - CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt - INTEGER :: i, j, l, ji, iml, jml - REAL :: phystep - REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp - REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder - REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf - REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil - -!--- Local variables for sea-ice reading: - LOGICAL :: read_mask - INTEGER :: iml_lic, jml_lic, isst(klon-2) - INTEGER :: fid, llm_tmp, ttm_tmp, itaul(1) - REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:) - REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:) - REAL :: date, lev(1), dummy - REAL :: flic_tmp(SIZE(masque,1),SIZE(masque,2)) - -!--- Arguments for conf_phys - LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats - REAL :: solarlong0, seuil_inversion, fact_cldcon, facttemps - LOGICAL :: ok_newmicro - INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs - REAL :: ratqsbas, ratqshaut, tau_ratqs - LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple - INTEGER :: flag_aerosol - LOGICAL :: flag_aerosol_strat - LOGICAL :: new_aod - REAL :: bl95_b0, bl95_b1 - INTEGER :: read_climoz !--- Read ozone climatology - REAL :: alp_offset - - deg2rad= pi/180.0 - iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml") - jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml") - -! Grid construction and miscellanous initializations. -!******************************************************************************* - CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & - callstats, & - solarlong0,seuil_inversion, & - fact_cldcon, facttemps,ok_newmicro,iflag_radia, & - iflag_cldcon, & - iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & - ok_ade, ok_aie, ok_cdnc, aerosol_couple, & - flag_aerosol, flag_aerosol_strat, new_aod, & - bl95_b0, bl95_b1, & - read_climoz, & - alp_offset) - - CALL phys_state_var_init(read_climoz) - - co2_ppm0 = co2_ppm !--- Initial atmospheric CO2 conc. from .def file - - rlat(1) = pi/2. - DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO - rlat(klon) = - pi/2. - rlat(:)=rlat(:)*(180.0/pi) - - rlon(1) = 0.0 - DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO - rlon(klon) = 0.0 - rlon(:)=rlon(:)*(180.0/pi) - -! Compute ground geopotential, sub-cells quantities and possibly the mask. -!******************************************************************************* - read_mask=ANY(masque/=-99999.); masque_tmp=masque - CALL start_init_orog(rlonv, rlatu, phis, masque_tmp) - WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt) - IF(.NOT.read_mask) THEN !--- Keep mask form orography - masque=masque_tmp - IF(prt_level>=1) THEN - WRITE(lunout,*)'BUILT MASK :' - WRITE(lunout,fmt) NINT(masque) - END IF - WHERE( masque(:,:)=1) & !--- Ozone climatology - CALL regr_lat_time_climoz(read_climoz) - -! Sub-surfaces initialization -!******************************************************************************* -!--- Read and interpolate on model T-grid soil fraction and soil ice fraction. - CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) - ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) - ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) - ALLOCATE( fraclic(iml_lic,jml_lic)) - CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & - & lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) - CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) - CALL flinclo(fid) - WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic - IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees - IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. - dlon_lic(:)=lon_lic(:,1) - dlat_lic(:)=lat_lic(1,:) - CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) ) - flic_tmp(iml,:)=flic_tmp(1,:) - -!--- To the physical grid - pctsrf(:,:) = 0. - CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic)) - -!--- Adequation with soil/sea mask - WHERE(pctsrf(:,is_lic)EPSFRA) THEN - IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN - pctsrf(ji,is_lic)=zmasq(ji) - pctsrf(ji,is_ter)=0. - ELSE - pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) - IF(pctsrf(ji,is_ter)0.) isst=1 - -!--- It is checked that the sub-surfaces sum is equal to 1. - ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) - IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points' - -! Write physical initial state -!******************************************************************************* - WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad - phystep = dtvr * FLOAT(iphysiq) - radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) ) - WRITE(lunout,*)'phystep =', phystep, radpas - -! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0 -!******************************************************************************* - DO i=1,nbsrf; ftsol(:,i) = tsol; END DO - DO i=1,nbsrf; snsrf(:,i) = sn; END DO - falb_dir(:,is_ter,:) = 0.08 - falb_dir(:,is_lic,:) = 0.6 - falb_dir(:,is_oce,:) = 0.5 - falb_dir(:,is_sic,:) = 0.6 - fevap(:,:) = 0. - DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO - DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO - rain_fall = 0. - snow_fall = 0. - solsw = 165. - sollw = -53. - t_ancien = 273.15 - q_ancien = 0. - agesno = 0. - - z0m(:,is_oce) = rugmer(:) - z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) - z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0) - z0m(:,is_sic) = 0.001 - z0h(:,:)=z0m(:,:) - - fder = 0.0 - clwcon = 0.0 - rnebcon = 0.0 - ratqs = 0.0 - run_off_lic_0 = 0.0 - rugoro = 0.0 - -! Before phyredem calling, surface modules and values to be saved in startphy.nc -! are initialized -!******************************************************************************* - dummy = 1.0 - pbl_tke(:,:,:) = 1.e-8 - zmax0(:) = 40. - f0(:) = 1.e-5 - sig1(:,:) = 0. - w01(:,:) = 0. - wake_deltat(:,:) = 0. - wake_deltaq(:,:) = 0. - wake_s(:) = 0. - wake_cstar(:) = 0. - wake_fip(:) = 0. - wake_pe = 0. - fm_therm = 0. - entr_therm = 0. - detr_therm = 0. - - CALL fonte_neige_init(run_off_lic_0) - CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil ) - CALL phyredem( "startphy.nc" ) - -! WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0' -! WRITE(lunout,*)'entree histclo' - CALL histclo() - -#endif -!#endif of #ifdef CPP_EARTH - -END SUBROUTINE etat0phys_netcdf -! -!------------------------------------------------------------------------------- - - -#ifdef CPP_EARTH -!------------------------------------------------------------------------------- -! -SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque) -! -!=============================================================================== -! Comment: -! This routine launch grid_noro, which computes parameters for SSO scheme as -! described in LOTT & MILLER (1997) and LOTT(1999). -!=============================================================================== - USE conf_dat_m, ONLY: conf_dat2d -! USE grid_atob_m, ONLY: rugsoro - USE grid_noro_m, ONLY: grid_noro - IMPLICIT NONE -!------------------------------------------------------------------------------- -! Arguments: - REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) - REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) -!------------------------------------------------------------------------------- -! Local variables: - CHARACTER(LEN=256) :: modname="start_init_orog" - CHARACTER(LEN=256) :: title="RELIEF" - INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) - REAL :: lev(1), date, dt - REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) - REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var (:,:) - REAL, ALLOCATABLE :: zmea0(:,:), zstd0(:,:), zsig0(:,:) - REAL, ALLOCATABLE :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:) -!------------------------------------------------------------------------------- - iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") - jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") - -!--- HIGH RESOLUTION OROGRAPHY - CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) - - ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) - CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& - lev, ttm_tmp, itau, date, dt, fid) - ALLOCATE(relief_hi(iml_rel,jml_rel)) - CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) - CALL flinclo(fid) - -!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS - ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) - lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad - lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad - -!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS - ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) - CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) - DEALLOCATE(lon_ini,lat_ini) - -!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro - WRITE(lunout,*) - WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***' - -!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES - ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation - ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy - ALLOCATE(zthe0(iml,jml)) !--- Highest slope orientation - ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights - -!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS - CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0, & - zsig0,zgam0,zthe0,zpic0,zval0,masque) - phis = phis * 9.81 - phis(iml,:) = phis(1,:) - -!--- COMPUTE SURFACE ROUGHNESS -! WRITE(lunout,*) -! WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***' -! ALLOCATE(tmp_var(iml-1,jml)) -! CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var) -! ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:) -! DEALLOCATE(tmp_var) - DEALLOCATE(relief_hi,lon_rad,lat_rad) - -!--- PUT QUANTITIES TO PHYSICAL GRID - CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0) - CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0) - CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0) - CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0) - CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0) - CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0) - CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0) - - -END SUBROUTINE start_init_orog -! -!------------------------------------------------------------------------------- - - -!------------------------------------------------------------------------------- -! -SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis) -! -!=============================================================================== -! Purpose: Compute tsol and qsol, knowing phis. -!=============================================================================== - IMPLICIT NONE -!------------------------------------------------------------------------------- -! Arguments: - REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) - REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) - LOGICAL, INTENT(IN) :: ibar - REAL, INTENT(IN) :: phis(:,:) ! dim (iml,jml) -!------------------------------------------------------------------------------- -! Local variables: - CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc" - REAL :: date, dt - INTEGER :: iml, jml, jml2, itau(1) - REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) - REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) - REAL, ALLOCATABLE :: ts(:,:), qs(:,:) -!------------------------------------------------------------------------------- - iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml") - jml=assert_eq(SIZE(lat_in),SIZE(phis,2), TRIM(modname)//" jml") - jml2=SIZE(lat_in2) - - WRITE(lunout,*)'Opening the surface analysis' - CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys) - WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys - - ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys)) - ALLOCATE(levphys_ini(llm_phys)) - CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys, & - lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys) - -!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS - ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys)) - lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad - lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad - - ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) - CALL get_var_phys('ST' ,ts) !--- SURFACE TEMPERATURE - CALL get_var_phys('CDSW',qs) !--- SOIL MOISTURE - CALL flinclo(fid_phys) - DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) - -!--- TSOL AND QSOL ON PHYSICAL GRID - ALLOCATE(tsol(klon)) - CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol) - CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol) - DEALLOCATE(ts,qs) - -CONTAINS - -!------------------------------------------------------------------------------- -! -SUBROUTINE get_var_phys(title,field) -! -!------------------------------------------------------------------------------- - USE conf_dat_m, ONLY: conf_dat2d - IMPLICIT NONE -!------------------------------------------------------------------------------- -! Arguments: - CHARACTER(LEN=*), INTENT(IN) :: title - REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) -!------------------------------------------------------------------------------- -! Local variables: - INTEGER :: tllm -!------------------------------------------------------------------------------- - SELECT CASE(title) - CASE('SP'); tllm=0 - CASE('ST','CDSW'); tllm=llm_phys - END SELECT - IF(ALLOCATED(field)) RETURN - ALLOCATE(field(iml,jml)); field(:,:)=0. - CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana) - CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar) - CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana, & - lon_in, lat_in, lon_in2, lat_in2, field) - -END SUBROUTINE get_var_phys -! -!------------------------------------------------------------------------------- -! -END SUBROUTINE start_init_phys -! -!------------------------------------------------------------------------------- - - -!------------------------------------------------------------------------------- -! -SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) -! -!------------------------------------------------------------------------------- - USE inter_barxy_m, ONLY: inter_barxy - USE grid_atob_m, ONLY: grille_m - IMPLICIT NONE -!------------------------------------------------------------------------------- -! Arguments: - CHARACTER(LEN=*), INTENT(IN) :: nam - LOGICAL, INTENT(IN) :: ibar, ibeg - REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) - REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) - REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) - REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) - REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) -!------------------------------------------------------------------------------- -! Local variables: - CHARACTER(LEN=256) :: modname="interp_startvar" - INTEGER :: ii, jj, i1, j1, j2 - REAL, ALLOCATABLE :: vtmp(:,:) -!------------------------------------------------------------------------------- - ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") - jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") - i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") - j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") - j2=SIZE(lat2) - ALLOCATE(vtmp(i1-1,j1)) - IF(ibar) THEN - IF(ibeg.AND.prt_level>1) THEN - WRITE(lunout,*)"--------------------------------------------------------" - WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" - WRITE(lunout,*)"--------------------------------------------------------" - END IF - CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) - ELSE - CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) - END IF - CALL gr_int_dyn(vtmp, varo, i1-1, j1) - -END SUBROUTINE interp_startvar -! -!------------------------------------------------------------------------------- - -#endif -!#endif of #ifdef CPP_EARTH - -END MODULE etat0phys -! -!******************************************************************************* -