Changeset 2336 for LMDZ5/trunk/libf/dynlonlat_phylonlat
- Timestamp:
- Jul 31, 2015, 7:22:21 PM (9 years ago)
- Location:
- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
r2331 r2336 1 1 PROGRAM ce0l 2 ! 3 ! Purpose: Calls etat0, creates initial states and limit_netcdf 4 ! 5 ! interbar=.T. for barycentric interpolation inter_barxy 6 ! extrap =.T. for data extrapolation, like for the SSTs when file does not 7 ! contain ocean points only. 8 ! oldice =.T. for old-style ice, obtained using grille_m (grid_atob). 9 ! masque is created in etat0, passed to limit to ensure consistancy. 10 11 USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS 12 USE etat0dyn, only: etat0dyn_netcdf 13 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 2 ! 3 !------------------------------------------------------------------------------- 4 ! Purpose: Initial states and boundary conditions files creation: 5 ! * start.nc for dynamics (using etat0dyn routine) 6 ! * startphy.nc for physics (using etat0phys routine) 7 ! * limit.nc for forced runs (using limit_netcdf routine) 8 !------------------------------------------------------------------------------- 9 ! Notes: 10 ! * extrap=.T. (default) for data extrapolation, like for the SSTs when file 11 ! does contain ocean points only. 12 ! * "masque" can be: 13 ! - read from file "o2a.nc" (for coupled runs). 14 ! - created in etat0phys or etat0dyn (for forced runs). 15 ! It is then passed to limit_netcdf to ensure consistancy. 16 !------------------------------------------------------------------------------- 14 17 USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo 15 16 USE etat0phys, only: etat0phys_netcdf 17 USE dimphy, only: KLON 18 USE infotrac, only: TYPE_TRAC, infotrac_init 18 USE control_mod, ONLY: day_step, dayref, nsplit_phys 19 USE etat0dyn, ONLY: etat0dyn_netcdf 20 USE etat0phys, ONLY: etat0phys_netcdf 21 USE limit, ONLY: limit_netcdf 22 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR 23 USE infotrac, ONLY: type_trac, infotrac_init 24 USE dimphy, ONLY: klon 19 25 USE test_disvert_m, ONLY: test_disvert 26 USE filtreg_mod, ONLY: inifilr 27 #ifdef inca 28 USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic 29 #endif 30 #ifdef CPP_PARA 31 USE mod_const_mpi, ONLY: init_const_mpi 32 USE parallel_lmdz, ONLY: init_parallel, mpi_rank, omp_rank, mpi_size 33 USE bands, ONLY: read_distrib, distrib_phys 34 USE mod_hallo, ONLY: init_mod_hallo 35 USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys 36 #endif 20 37 21 38 IMPLICIT NONE 22 39 23 ! Local variables: 40 !------------------------------------------------------------------------------- 41 ! Local variables: 24 42 include "dimensions.h" 25 43 include "paramet.h" 26 include "comgeom .h"44 include "comgeom2.h" 27 45 include "comconst.h" 28 46 include "comvert.h" … … 30 48 include "temps.h" 31 49 include "logic.h" 32 REAL :: masque(iip1, 33 REAL :: phis (iip1, 50 REAL :: masque(iip1,jjp1) !--- CONTINENTAL MASK 51 REAL :: phis (iip1,jjp1) !--- GROUND GEOPOTENTIAL 34 52 CHARACTER(LEN=256) :: modname, fmt, calnd !--- CALENDAR TYPE 35 53 LOGICAL :: use_filtre_fft 36 LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.37 38 54 LOGICAL, PARAMETER :: extrap=.FALSE. 55 56 !--- Local variables for ocean mask reading: 39 57 INTEGER :: nid_o2a, iml_omask, jml_omask, j 40 58 INTEGER :: fid, iret, llm_tmp, ttm_tmp, itaul(1) 41 REAL, ALLOCATABLE :: lon_omask(:, :), dlon_omask(:), ocemask(:,:)42 REAL, ALLOCATABLE :: lat_omask(:, :), dlat_omask(:), ocetmp (:,:)59 REAL, ALLOCATABLE :: lon_omask(:,:), dlon_omask(:), ocemask(:,:) 60 REAL, ALLOCATABLE :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:) 43 61 REAL :: date, lev(1) 44 45 !---------------------------------------------------------------------- 62 !------------------------------------------------------------------------------- 46 63 modname="ce0l" 47 64 48 65 !--- Constants 49 66 pi = 4. * ATAN(1.) 50 67 rad = 6371229. … … 59 76 60 77 CALL conf_gcm( 99, .TRUE. ) 61 62 78 dtvr = daysec/REAL(day_step) 63 WRITE(lunout, *)'dtvr', dtvr 64 79 WRITE(lunout,*)'dtvr',dtvr 65 80 CALL iniconst() 66 81 CALL inigeom() 67 82 83 !--- Calendar choice 68 84 #ifdef CPP_IOIPSL 69 85 calnd='gregorian' 70 86 SELECT CASE(calend) 71 CASE('earth_360d') 72 CALL ioconf_calendar('360d') 73 calnd='with 360 days/year' 74 CASE('earth_365d') 75 CALL ioconf_calendar('noleap') 76 calnd='with no leap year' 77 CASE('earth_366d') 78 CALL ioconf_calendar('366d') 79 calnd='with leap years only' 80 CASE('gregorian') 81 CALL ioconf_calendar('gregorian') 82 CASE('standard') 83 CALL ioconf_calendar('gregorian') 84 CASE('julian') 85 CALL ioconf_calendar('julian') 86 calnd='julian' 87 CASE('proleptic_gregorian') 88 CALL ioconf_calendar('gregorian') 89 !--- DC Bof... => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian 90 CASE DEFAULT 91 CALL abort_gcm('ce0l', 'Bad choice for calendar', 1) 87 CASE('earth_360d');CALL ioconf_calendar('360d'); calnd='with 360 days/year' 88 CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year' 89 CASE('earth_366d');CALL ioconf_calendar('366d'); calnd='with leap years only' 90 CASE('gregorian'); CALL ioconf_calendar('gregorian') 91 CASE('standard'); CALL ioconf_calendar('gregorian') 92 CASE('julian'); CALL ioconf_calendar('julian'); calnd='julian' 93 CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian') 94 !--- DC Bof... => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian 95 CASE DEFAULT 96 CALL abort_gcm('ce0l','Bad choice for calendar',1) 92 97 END SELECT 93 WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd) 94 #endif 95 96 use_filtre_fft=.FALSE. 97 CALL getin('use_filtre_fft', use_filtre_fft) 98 IF(use_filtre_fft) THEN 99 WRITE(lunout, *)"FFT filter not available for sequential dynamics." 100 WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used." 101 ENDIF 102 103 !--- LAND MASK. TWO CASES: 104 ! 1) read from ocean model file "o2a.nc" (coupled runs) 105 ! 2) computed from topography file "Relief.nc" (masque(:, :)=-99999.) 106 ! Coupled simulations (case 1) use the ocean model mask to compute the 107 ! weights to ensure ocean fractions are the same for atmosphere and ocean. 108 109 IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN 110 WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found' 111 WRITE(lunout, *)'Forced run.' 112 masque(:, :)=-99999. 113 ELSE 114 iret=NF90_CLOSE(nid_o2a) 115 WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found' 116 WRITE(lunout, *)'Coupled run.' 117 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 118 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 119 WRITE(lunout, *)'Mismatching dimensions for ocean mask' 120 WRITE(lunout, *)'iim = ', iim , ' iml_omask = ', iml_omask 121 WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask 122 CALL abort_gcm(modname, '', 1) 123 END IF 124 ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim )) 125 ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1)) 126 CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, & 127 lat_omask, lev, ttm_tmp, itaul, date, dt, fid) 128 CALL flinget(fid, "OceMask", iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp) 129 CALL flinclo(fid) 130 dlon_omask(1:iim ) = lon_omask(1:iim, 1) 131 dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1) 132 ocemask = ocetmp 133 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 134 DO j=1, jjp1 135 ocemask(:, j) = ocetmp(:, jjp1-j+1) 136 END DO 137 END IF 138 DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask) 139 IF(prt_level>=1) THEN 140 WRITE(fmt, "(i4, 'i1)')")iim 141 fmt='('//ADJUSTL(fmt) 142 WRITE(lunout, *)'OCEAN MASK :' 143 WRITE(lunout, fmt) NINT(ocemask) 144 END IF 145 masque(1:iim, :)=1.-ocemask(:, :) 146 masque(iip1 , :)=masque(1, :) 147 DEALLOCATE(ocemask) 148 END IF 149 phis(:, :)=-99999. 150 151 CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/)) 152 WRITE(lunout, *)'---> klon=', klon 153 154 call infotrac_init 155 CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, & 156 rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys) 157 158 IF(pressure_exner) CALL test_disvert 159 98 WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd) 99 #endif 100 101 #ifdef CPP_PARA 102 !--- Physical grid + parallel initializations 103 CALL init_const_mpi() 104 CALL init_parallel() 105 CALL read_distrib() 106 CALL init_mod_hallo() 107 CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys) 108 CALL init_interface_dyn_phys() 109 #else 110 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 111 #endif 112 WRITE(lunout,*)'---> klon=',klon 113 114 !--- Tracers initializations 160 115 IF (type_trac == 'inca') THEN 161 116 #ifdef INCA 162 CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday) 163 CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0) 164 WRITE(lunout, *)'nbtr =' , nbtr 165 #endif 166 END IF 117 CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,& 118 nbsrf,is_oce,is_sic,is_ter,is_lic,calend) 119 CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,& 120 COMM_LMDZ) 121 WRITE(lunout,*)'nbtr =' , nbtr 122 #endif 123 END IF 124 CALL infotrac_init() 125 126 CALL inifilr() 127 CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, & 128 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) 129 IF(pressure_exner) CALL test_disvert 130 131 #ifdef CPP_PARA 132 IF (mpi_rank==0.AND.omp_rank==0) THEN 133 #endif 134 use_filtre_fft=.FALSE. 135 CALL getin('use_filtre_fft',use_filtre_fft) 136 IF(use_filtre_fft) THEN 137 WRITE(lunout,*)"FFT filter not available for sequential dynamics." 138 WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used." 139 ENDIF 140 141 !--- LAND MASK. TWO CASES: 142 ! 1) read from ocean model file "o2a.nc" (coupled runs) 143 ! 2) computed from topography file "Relief.nc" (masque(:,:)=-99999.) 144 ! Coupled simulations (case 1) use the ocean model mask to compute the 145 ! weights to ensure ocean fractions are the same for atmosphere and ocean. 146 !******************************************************************************* 147 IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN 148 WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found' 149 WRITE(lunout,*)'Forced run.' 150 masque(:,:)=-99999. 151 ELSE 152 iret=NF90_CLOSE(nid_o2a) 153 WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found' 154 WRITE(lunout,*)'Coupled run.' 155 CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a) 156 IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN 157 WRITE(lunout,*)'Mismatching dimensions for ocean mask' 158 WRITE(lunout,*)'iim = ',iim ,' iml_omask = ',iml_omask 159 WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask 160 CALL abort_gcm(modname,'',1) 161 END IF 162 ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim )) 163 ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1)) 164 CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev, & 165 & ttm_tmp,itaul,date,dt,fid) 166 CALL flinget(fid, "OceMask", iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp) 167 CALL flinclo(fid) 168 dlon_omask(1:iim ) = lon_omask(1:iim,1) 169 dlat_omask(1:jjp1) = lat_omask(1,1:jjp1) 170 ocemask = ocetmp 171 IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN 172 DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO 173 END IF 174 DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask) 175 IF(prt_level>=1) THEN 176 WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt) 177 WRITE(lunout,*)'OCEAN MASK :' 178 WRITE(lunout,fmt) NINT(ocemask) 179 END IF 180 masque(1:iim,:)=1.-ocemask(:,:) 181 masque(iip1 ,:)=masque(1,:) 182 DEALLOCATE(ocemask) 183 END IF 184 phis(:,:)=-99999. 185 167 186 IF(ok_etat0) THEN 168 WRITE(lunout, '(//)') 169 WRITE(lunout, *) ' ************************ ' 170 WRITE(lunout, *) ' *** etat0phy_netcdf *** ' 171 WRITE(lunout, *) ' ************************ ' 172 WRITE(lunout, '(//)') 173 WRITE(lunout, *) ' interbar = ', interbar 174 CALL etat0phys_netcdf(interbar, masque, phis) 175 END IF 176 177 IF(ok_etat0) THEN 178 WRITE(lunout, '(//)') 179 WRITE(lunout, *) ' ************************ ' 180 WRITE(lunout, *) ' *** etat0dyn_netcdf *** ' 181 WRITE(lunout, *) ' ************************ ' 182 WRITE(lunout, '(//)') 183 WRITE(lunout, *) ' interbar = ', interbar 184 CALL etat0dyn_netcdf(interbar, masque, phis) 187 WRITE(lunout,'(//)') 188 WRITE(lunout,*) ' ************************ ' 189 WRITE(lunout,*) ' *** etat0phy_netcdf *** ' 190 WRITE(lunout,*) ' ************************ ' 191 CALL etat0phys_netcdf(masque,phis) 192 WRITE(lunout,'(//)') 193 WRITE(lunout,*) ' ************************ ' 194 WRITE(lunout,*) ' *** etat0dyn_netcdf *** ' 195 WRITE(lunout,*) ' ************************ ' 196 CALL etat0dyn_netcdf(masque,phis) 185 197 END IF 186 198 187 199 IF(ok_limit) THEN 188 WRITE(lunout, '(//)') 189 WRITE(lunout, *) ' ********************* ' 190 WRITE(lunout, *) ' *** Limit_netcdf *** ' 191 WRITE(lunout, *) ' ********************* ' 192 WRITE(lunout, '(//)') 193 CALL limit_netcdf(interbar, extrap, oldice, masque) 194 END IF 195 196 WRITE(lunout, '(//)') 197 WRITE(lunout, *) ' *************************** ' 198 WRITE(lunout, *) ' *** grilles_gcm_netcdf *** ' 199 WRITE(lunout, *) ' *************************** ' 200 WRITE(lunout, '(//)') 201 CALL grilles_gcm_netcdf_sub(masque, phis) 200 WRITE(lunout,'(//)') 201 WRITE(lunout,*) ' ********************* ' 202 WRITE(lunout,*) ' *** Limit_netcdf *** ' 203 WRITE(lunout,*) ' ********************* ' 204 WRITE(lunout,'(//)') 205 CALL limit_netcdf(masque,phis,extrap) 206 END IF 207 208 WRITE(lunout,'(//)') 209 WRITE(lunout,*) ' *************************** ' 210 WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' 211 WRITE(lunout,*) ' *************************** ' 212 WRITE(lunout,'(//)') 213 CALL grilles_gcm_netcdf_sub(masque,phis) 214 215 #ifdef CPP_PARA 216 END IF 217 #endif 202 218 203 219 END PROGRAM ce0l 220 ! 221 !------------------------------------------------------------------------------- -
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
r2334 r2336 12 12 ! routine (to be called after restget): 13 13 ! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,& 14 ! champ, lon_in2, lat_in2 , ibar)14 ! champ, lon_in2, lat_in2) 15 15 ! 16 16 ! * Variables should have the following names in the NetCDF files: … … 36 36 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo 37 37 USE assert_eq_m, ONLY: assert_eq 38 #ifdef CPP_PHYS39 38 USE indice_sol_mod, ONLY: epsfra 40 #endif41 39 IMPLICIT NONE 42 40 … … 54 52 include "serre.h" 55 53 REAL, SAVE :: deg2rad 56 #ifndef CPP_PHYS57 REAL, SAVE :: epsfra= 1.E-558 #endif59 54 INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn 60 55 REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:) 61 56 CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc' 62 CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'63 57 64 58 CONTAINS … … 66 60 !------------------------------------------------------------------------------- 67 61 ! 68 SUBROUTINE etat0dyn_netcdf( ib,masque, phis)62 SUBROUTINE etat0dyn_netcdf(masque, phis) 69 63 ! 70 64 !------------------------------------------------------------------------------- … … 73 67 ! Notes: 1) This routine is designed to work for Earth 74 68 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 75 ! Otherwise: read it from ocean model file (coupled run) orcompute it.69 ! Otherwise: compute it. 76 70 !------------------------------------------------------------------------------- 77 71 USE control_mod 78 #ifdef CPP_PHYS79 72 USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz 80 73 USE regr_pr_o3_m, ONLY: regr_pr_o3 81 74 USE press_coefoz_m, ONLY: press_coefoz 82 #endif83 75 USE exner_hyb_m, ONLY: exner_hyb 84 76 USE exner_milieu_m, ONLY: exner_milieu 85 USE infotrac, only: NQTOT, TNAME77 USE infotrac, ONLY: nqtot, tname 86 78 USE filtreg_mod 87 79 IMPLICIT NONE 88 80 !------------------------------------------------------------------------------- 89 81 ! Arguments: 90 LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation91 82 REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask 92 83 REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential … … 111 102 deg2rad = pi/180.0 112 103 113 ! Initializations for tracers and filter114 !*******************************************************************************115 CALL inifilr()116 117 104 ! Compute ground geopotential and possibly the mask. 118 105 !******************************************************************************* 119 106 masque_tmp(:,:)=masque(:,:) 120 CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)121 107 WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt) 122 108 IF(ALL(masque==-99999.)) THEN !--- KEEP NEW MASK … … 132 118 ! Compute psol AND tsol, knowing phis. 133 119 !******************************************************************************* 134 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib,phis, psol)120 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) 135 121 136 122 ! Mid-levels pressure computation … … 147 133 !******************************************************************************* 148 134 uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0 149 CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv ,ib)135 CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv) 150 136 CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, & 151 & rlonu,rlatu(:jjm) ,ib)152 CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv ,ib)137 & rlonu,rlatu(:jjm)) 138 CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv) 153 139 tpot(:,:,:)=t3d(:,:,:) 154 CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv ,ib)140 CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv) 155 141 156 142 WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:)) … … 165 151 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 166 152 qd (:,:,:) = 0.0 167 CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv ,ib)153 CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv) 168 154 ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:) 169 155 CALL flinclo(fid_dyn) 170 156 171 #ifdef CPP_PHYS172 157 ! Parameterization of ozone chemistry: 173 158 !******************************************************************************* … … 180 165 q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29. !--- Mole->mass fraction 181 166 END IF 182 183 #endif184 167 q3d(iip1,:,:,:)=q3d(1,:,:,:) 185 186 ! Intermediate computation187 !*******************************************************************************188 CALL massdair(p3d,masse)189 WRITE(lunout,*)' ALPHAX ',alphax190 DO l=1,llm191 xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l)192 xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)193 xpn=SUM(xppn)/apoln194 xps=SUM(xpps)/apols195 masse(:,1 ,l)=xpn196 masse(:,jjp1,l)=xps197 END DO198 168 199 169 ! Writing … … 215 185 CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) 216 186 WRITE(lunout,*)'sortie geopot' 217 CALL caldyn0( itau, uvent, vvent, tpot, psol, pk, phis, &218 phi, w, time+iday-dayref)187 CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & 188 phi, w, pbaru, pbarv, time+iday-dayref) 219 189 WRITE(lunout,*)'sortie caldyn0' 190 #ifdef CPP_PARA 191 CALL dynredem0_loc( "start.nc", dayref, phis) 192 #else 220 193 CALL dynredem0( "start.nc", dayref, phis) 194 #endif 221 195 WRITE(lunout,*)'sortie dynredem0' 196 #ifdef CPP_PARA 197 CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 198 #else 222 199 CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 200 #endif 223 201 WRITE(lunout,*)'sortie dynredem1' 224 202 CALL histclo() … … 232 210 ! 233 211 SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& 234 champ, lon_in2, lat_in2 , ibar)212 champ, lon_in2, lat_in2) 235 213 !------------------------------------------------------------------------------- 236 214 IMPLICIT NONE … … 252 230 REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) 253 231 REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) 254 LOGICAL, INTENT(IN) :: ibar255 232 !------------------------------------------------------------------------------- 256 233 ! Local variables: … … 287 264 !--- INTERPOLATE 3D FIELD IF NEEDED 288 265 IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & 289 lat_in2,pls,champ ,ibar)266 lat_in2,pls,champ) 290 267 291 268 !--- COMPUTE THE REQUIRED FILED … … 317 294 !------------------------------------------------------------------------------- 318 295 ! 319 SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque) 320 ! 321 !------------------------------------------------------------------------------- 322 USE conf_dat_m, ONLY: conf_dat2d 323 IMPLICIT NONE 324 !=============================================================================== 325 ! Purpose: Compute "phis" just like it would be in start_init_orog. 326 !=============================================================================== 327 ! Arguments: 328 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) 329 REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) 330 !------------------------------------------------------------------------------- 331 ! Local variables: 332 CHARACTER(LEN=256) :: modname="start_init_orog0" 333 CHARACTER(LEN=256) :: title="RELIEF" 334 INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) 335 REAL :: lev(1), date, dt 336 REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) 337 REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:) 338 !------------------------------------------------------------------------------- 339 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") 340 jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") 341 IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1) 342 IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1) 343 pi=2.0*ASIN(1.0); deg2rad=pi/180.0 344 IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN 345 346 !--- HIGH RESOLUTION OROGRAPHY 347 CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) 348 349 ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) 350 CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& 351 lev, ttm_tmp, itau, date, dt, fid) 352 ALLOCATE(relief_hi(iml_rel,jml_rel)) 353 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) 354 CALL flinclo(fid) 355 356 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 357 ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) 358 lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad 359 lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad 360 361 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 362 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 363 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) 364 DEALLOCATE(lon_ini,lat_ini) 365 366 !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0 367 WRITE(lunout,*) 368 WRITE(lunout,*)'*** Compute surface geopotential ***' 369 370 !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS 371 CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque) 372 phis = phis * 9.81 373 phis(iml,:) = phis(1,:) 374 DEALLOCATE(relief_hi,lon_rad,lat_rad) 375 376 END SUBROUTINE start_init_orog0 377 ! 378 !------------------------------------------------------------------------------- 379 380 381 !------------------------------------------------------------------------------- 382 ! 383 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask) 384 ! 385 !=============================================================================== 386 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics 387 ! without any call to physics subroutines. 388 !=============================================================================== 389 IMPLICIT NONE 390 !------------------------------------------------------------------------------- 391 ! Arguments: 392 REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) 393 REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp) 394 REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) 395 REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) 396 REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar) 397 !------------------------------------------------------------------------------- 398 ! Local variables: 399 CHARACTER(LEN=256) :: modname="grid_noro0" 400 REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) 401 REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2) 402 REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar) 403 REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar) 404 REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax) 405 REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax) 406 REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax) 407 LOGICAL :: masque_lu 408 INTEGER :: i, ii, imdp, imar, iext 409 INTEGER :: j, jj, jmdp, jmar, nn 410 REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest 411 REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue 412 !------------------------------------------------------------------------------- 413 imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") 414 jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") 415 imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1 416 jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar") 417 IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1) 418 IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1) 419 iext=imdp/10 420 xpi = ACOS(-1.) 421 rad = 6371229. 422 423 !--- ARE WE USING A READ MASK ? 424 masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 425 WRITE(lunout,*)'Masque lu: ',masque_lu 426 427 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: 428 ALLOCATE(xusn(imdp+2*iext)) 429 xusn(1 +iext:imdp +iext)=xd(:) 430 xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi 431 xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi 432 433 ALLOCATE(yusn(jmdp+2)) 434 yusn(1 )=yd(1) +(yd(1) -yd(2)) 435 yusn(2:jmdp+1)=yd(:) 436 yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) 437 438 ALLOCATE(zusn(imdp+2*iext,jmdp+2)) 439 zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) 440 zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) 441 zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) 442 zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) 443 zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) 444 zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) 445 zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) 446 447 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) 448 ALLOCATE(a(imar+1),b(imar+1)) 449 b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 450 b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 451 a(1)=x(1)-(x(2)-x(1))/2.0 452 a(2:imar+1)= b(1:imar) 453 454 ALLOCATE(c(jmar),d(jmar)) 455 d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 456 d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 457 c(1)=y(1)-(y(2)-y(1))/2.0 458 c(2:jmar)=d(1:jmar-1) 459 460 !--- INITIALIZATIONS: 461 ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0 462 ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0 463 464 !--- SUMMATION OVER GRIDPOINT AREA 465 zleny=xpi/REAL(jmdp)*rad 466 xincr=xpi/REAL(jmdp)/2. 467 ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. 468 ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. 469 DO ii = 1, imar+1 470 DO jj = 1, jmar 471 DO j = 2,jmdp+1 472 zlenx =zleny *COS(yusn(j)) 473 zbordnor=(xincr+c(jj)-yusn(j))*rad 474 zbordsud=(xincr-d(jj)+yusn(j))*rad 475 weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny)) 476 IF(weighy/=0) THEN 477 DO i = 2, imdp+2*iext-1 478 zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j)) 479 zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j)) 480 weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx)) 481 IF(weighx/=0)THEN 482 num_tot(ii,jj)=num_tot(ii,jj)+1.0 483 IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 484 weight(ii,jj)=weight(ii,jj)+weighx*weighy 485 zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN 486 END IF 487 END DO 488 END IF 489 END DO 490 END DO 491 END DO 492 493 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME 494 IF(.NOT.masque_lu) THEN 495 WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:) 496 END IF 497 nn=COUNT(weight(:,1:jmar-1)==0.0) 498 IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn 499 WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:) 500 501 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS) 502 ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0 503 WHERE(mask>=0.1) mask_tmp = 1. 504 WHERE(weight(:,:)/=0.0) 505 zphi(:,:)=mask_tmp(:,:)*zmea(:,:) 506 zmea(:,:)=mask_tmp(:,:)*zmea(:,:) 507 END WHERE 508 WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) 509 510 !--- Values at poles 511 zphi(imar+1,:)=zphi(1,:) 512 513 zweinor=SUM(weight(1:imar, 1),DIM=1) 514 zweisud=SUM(weight(1:imar,jmar),DIM=1) 515 zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1) 516 zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1) 517 zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud 518 519 END SUBROUTINE grid_noro0 520 ! 521 !------------------------------------------------------------------------------- 522 523 524 !------------------------------------------------------------------------------- 525 ! 526 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol) 296 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol) 527 297 ! 528 298 !------------------------------------------------------------------------------- … … 534 304 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 535 305 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 536 LOGICAL, INTENT(IN) :: ibar537 306 REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) 538 307 REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) … … 606 375 ALLOCATE(field(iml,jml)) 607 376 CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) 608 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)609 CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,&610 lon_in, lat_in, lon_in2, lat_in2, field)377 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 378 CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, & 379 lon_in, lat_in, lon_in2, lat_in2, field) 611 380 ELSE IF(SIZE(field)/=SIZE(z)) THEN 612 381 msg='The '//TRIM(msg)//' field we have does not have the right size' … … 625 394 !------------------------------------------------------------------------------- 626 395 ! 627 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d ,ibar)396 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d) 628 397 ! 629 398 !------------------------------------------------------------------------------- … … 639 408 REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) 640 409 REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) 641 LOGICAL, INTENT(IN) :: ibar642 410 !------------------------------------------------------------------------------- 643 411 ! Local variables: … … 667 435 ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) 668 436 CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & 669 lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)437 lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) 670 438 DEALLOCATE(lon_ini, lat_ini) 671 439 … … 673 441 ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) 674 442 DO il = 1,llm_dyn 675 CALL interp_startvar(var,i bar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),&443 CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & 676 444 lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) 677 445 END DO … … 708 476 !------------------------------------------------------------------------------- 709 477 ! 710 SUBROUTINE interp_startvar(nam,ib ar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)478 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) 711 479 ! 712 480 !------------------------------------------------------------------------------- 713 481 USE inter_barxy_m, ONLY: inter_barxy 714 USE grid_atob_m, ONLY: grille_m715 482 IMPLICIT NONE 716 483 !------------------------------------------------------------------------------- 717 484 ! Arguments: 718 485 CHARACTER(LEN=*), INTENT(IN) :: nam 719 LOGICAL, INTENT(IN) :: ib ar, ibeg486 LOGICAL, INTENT(IN) :: ibeg 720 487 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 721 488 REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) … … 735 502 j2=SIZE(lat2) 736 503 ALLOCATE(vtmp(i1-1,j1)) 737 IF(ibar) THEN 738 IF(ibeg.AND.prt_level>1) THEN 739 WRITE(lunout,*)"---------------------------------------------------------" 740 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 741 WRITE(lunout,*)"---------------------------------------------------------" 742 END IF 743 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 744 ELSE 745 CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) 746 END IF 504 IF(ibeg.AND.prt_level>1) THEN 505 WRITE(lunout,*)"---------------------------------------------------------" 506 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 507 WRITE(lunout,*)"---------------------------------------------------------" 508 END IF 509 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 747 510 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 748 511 -
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
r2320 r2336 36 36 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo 37 37 USE assert_eq_m, ONLY: assert_eq 38 USE dimphy 38 USE dimphy, ONLY: klon 39 USE conf_dat_m, ONLY: conf_dat2d 39 40 USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, & 40 41 rlon, solsw, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, & … … 51 52 include "paramet.h" 52 53 include "comgeom2.h" 53 include "comvert.h"54 54 include "comconst.h" 55 55 include "dimsoil.h" 56 56 include "temps.h" 57 include "comdissnew.h"58 include "serre.h"59 57 include "clesphys.h" 60 58 REAL, SAVE :: deg2rad 61 59 REAL, SAVE, ALLOCATABLE :: tsol(:) 62 ! REAL, SAVE, ALLOCATABLE :: rugo(:,:) ! ??? COMPUTED BUT NOT USED ???63 60 INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys 64 61 REAL, ALLOCATABLE, SAVE :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:) 65 CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc" 66 CHARACTER(LEN=256), PARAMETER :: title="RELIEF" 62 CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF" 63 CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc", psrfvar="SP" 64 CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW", tsrfvar="ST" 67 65 68 66 … … 72 70 !------------------------------------------------------------------------------- 73 71 ! 74 SUBROUTINE etat0phys_netcdf( ib,masque, phis)72 SUBROUTINE etat0phys_netcdf(masque, phis) 75 73 ! 76 74 !------------------------------------------------------------------------------- 77 75 ! Purpose: Creates initial states 78 76 !------------------------------------------------------------------------------- 79 ! Note: This routine is designed to work for Earth 77 ! Notes: 1) This routine is designed to work for Earth 78 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 79 ! Otherwise: compute it. 80 80 !------------------------------------------------------------------------------- 81 81 USE control_mod … … 84 84 USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz 85 85 USE indice_sol_mod 86 USE conf_phys_m, ONLY: conf_phys 87 USE exner_hyb_m, ONLY: exner_hyb 88 USE exner_milieu_m, ONLY: exner_milieu 89 USE test_disvert_m, ONLY: test_disvert 90 USE grid_atob_m, ONLY: grille_m 86 USE conf_phys_m, ONLY: conf_phys 87 USE init_ssrf_m, ONLY: start_init_subsurf 91 88 IMPLICIT NONE 92 89 !------------------------------------------------------------------------------- 93 90 ! Arguments: 94 LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation95 91 REAL, INTENT(INOUT) :: masque(:,:) !--- Land mask dim(iip1,jjp1) 96 92 REAL, INTENT(INOUT) :: phis (:,:) !--- Ground geopotential dim(iip1,jjp1) … … 99 95 CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt 100 96 INTEGER :: i, j, l, ji, iml, jml 101 REAL :: phystep 97 LOGICAL :: read_mask 98 REAL :: phystep, dummy 102 99 REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp 103 100 REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder 104 101 REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf 105 102 REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil 106 107 !--- Local variables for sea-ice reading:108 LOGICAL :: read_mask109 INTEGER :: iml_lic, jml_lic, isst(klon-2)110 INTEGER :: fid, llm_tmp, ttm_tmp, itaul(1)111 REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)112 REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)113 REAL :: date, lev(1), dummy114 REAL :: flic_tmp(SIZE(masque,1),SIZE(masque,2))115 103 116 104 !--- Arguments for conf_phys … … 132 120 jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml") 133 121 134 ! Grid construction and miscellanous initializations.122 ! Physics configuration 135 123 !******************************************************************************* 136 124 CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & … … 145 133 read_climoz, & 146 134 alp_offset) 147 148 135 CALL phys_state_var_init(read_climoz) 149 136 150 co2_ppm0 = co2_ppm !--- Initial atmospheric CO2 conc. from .def file 151 152 rlat(1) = pi/2. 153 DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j); END DO 154 rlat(klon) = - pi/2. 155 rlat(:)=rlat(:)*(180.0/pi) 156 157 rlon(1) = 0.0 158 DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO 159 rlon(klon) = 0.0 160 rlon(:)=rlon(:)*(180.0/pi) 137 !--- Initial atmospheric CO2 conc. from .def file 138 co2_ppm0 = co2_ppm 161 139 162 140 ! Compute ground geopotential, sub-cells quantities and possibly the mask. … … 174 152 WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1. 175 153 END IF 176 CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid154 CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid 177 155 178 156 ! Compute tsol and qsol on physical grid, knowing phis on 2D grid. 179 157 !******************************************************************************* 180 CALL start_init_phys(rlon v, rlatu, rlonu, rlatv, ib, phis)158 CALL start_init_phys(rlonu, rlatv, phis) 181 159 182 160 ! Some initializations. … … 188 166 CALL regr_lat_time_climoz(read_climoz) 189 167 190 ! Sub-surfaces initialization 191 !******************************************************************************* 192 !--- Read and interpolate on model T-grid soil fraction and soil ice fraction. 193 CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid) 194 ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic)) 195 ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic)) 196 ALLOCATE( fraclic(iml_lic,jml_lic)) 197 CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, & 198 & lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid) 199 CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic) 200 CALL flinclo(fid) 201 WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic 202 IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees 203 IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180. 204 dlon_lic(:)=lon_lic(:,1) 205 dlat_lic(:)=lat_lic(1,:) 206 CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) ) 207 flic_tmp(iml,:)=flic_tmp(1,:) 208 209 !--- To the physical grid 210 pctsrf(:,:) = 0. 211 CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic)) 212 213 !--- Adequation with soil/sea mask 214 WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 215 WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0. 216 pctsrf(:,is_ter)=zmasq(:) 217 DO ji=1,klon 218 IF(zmasq(ji)>EPSFRA) THEN 219 IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN 220 pctsrf(ji,is_lic)=zmasq(ji) 221 pctsrf(ji,is_ter)=0. 222 ELSE 223 pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic) 224 IF(pctsrf(ji,is_ter)<EPSFRA) THEN 225 pctsrf(ji,is_ter)=0. 226 pctsrf(ji,is_lic)=zmasq(ji) 227 END IF 228 END IF 229 END IF 230 END DO 231 232 !--- Sub-surface ocean and sea ice (sea ice set to zero for start). 233 pctsrf(:,is_oce)=(1.-zmasq(:)) 234 WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0. 235 IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:) 236 isst=0 237 WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1 238 239 !--- It is checked that the sub-surfaces sum is equal to 1. 240 ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA) 241 IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points' 168 ! Sub-surfaces initialization. 169 !******************************************************************************* 170 CALL start_init_subsurf(read_mask) 242 171 243 172 ! Write physical initial state … … 321 250 ! described in LOTT & MILLER (1997) and LOTT(1999). 322 251 !=============================================================================== 323 USE conf_dat_m, ONLY: conf_dat2d324 ! USE grid_atob_m, ONLY: rugsoro325 252 USE grid_noro_m, ONLY: grid_noro 326 253 IMPLICIT NONE … … 331 258 !------------------------------------------------------------------------------- 332 259 ! Local variables: 333 CHARACTER(LEN=256) :: modname="start_init_orog" 334 CHARACTER(LEN=256) :: title="RELIEF" 260 CHARACTER(LEN=256) :: modname 335 261 INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) 336 262 REAL :: lev(1), date, dt … … 340 266 REAL, ALLOCATABLE :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:) 341 267 !------------------------------------------------------------------------------- 268 modname="start_init_orog" 342 269 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") 343 270 jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") … … 350 277 lev, ttm_tmp, itau, date, dt, fid) 351 278 ALLOCATE(relief_hi(iml_rel,jml_rel)) 352 CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)279 CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi) 353 280 CALL flinclo(fid) 354 281 … … 360 287 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 361 288 ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) 362 CALL conf_dat2d( title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)289 CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.) 363 290 DEALLOCATE(lon_ini,lat_ini) 364 291 … … 378 305 phis = phis * 9.81 379 306 phis(iml,:) = phis(1,:) 380 381 !--- COMPUTE SURFACE ROUGHNESS382 ! WRITE(lunout,*)383 ! WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'384 ! ALLOCATE(tmp_var(iml-1,jml))385 ! CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)386 ! ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)387 ! DEALLOCATE(tmp_var)388 307 DEALLOCATE(relief_hi,lon_rad,lat_rad) 389 308 … … 405 324 !------------------------------------------------------------------------------- 406 325 ! 407 SUBROUTINE start_init_phys(lon_in,lat_in, lon_in2,lat_in2,ibar,phis)326 SUBROUTINE start_init_phys(lon_in,lat_in,phis) 408 327 ! 409 328 !=============================================================================== … … 413 332 !------------------------------------------------------------------------------- 414 333 ! Arguments: 415 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 416 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 417 LOGICAL, INTENT(IN) :: ibar 334 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml2) 418 335 REAL, INTENT(IN) :: phis(:,:) ! dim (iml,jml) 419 336 !------------------------------------------------------------------------------- 420 337 ! Local variables: 421 CHARACTER(LEN=256) :: modname ="start_init_phys", physfname="ECPHY.nc"338 CHARACTER(LEN=256) :: modname 422 339 REAL :: date, dt 423 340 INTEGER :: iml, jml, jml2, itau(1) … … 426 343 REAL, ALLOCATABLE :: ts(:,:), qs(:,:) 427 344 !------------------------------------------------------------------------------- 428 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")429 jml=assert_eq(SIZE(lat_in),SIZE(phis,2), TRIM(modname)//" jml")430 jml 2=SIZE(lat_in2)345 modname="start_init_phys" 346 iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml") 347 jml=SIZE(phis,2); jml2=SIZE(lat_in) 431 348 432 349 WRITE(lunout,*)'Opening the surface analysis' 433 CALL flininfo(phy sfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)434 WRITE(lunout,*) 'Values read: ', 350 CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys) 351 WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys 435 352 436 353 ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys)) 437 354 ALLOCATE(levphys_ini(llm_phys)) 438 CALL flinopen(phy sfname, .FALSE., iml_phys, jml_phys, llm_phys, &355 CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys, & 439 356 lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys) 440 357 … … 445 362 446 363 ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys)) 447 CALL get_var_phys( 'ST',ts) !--- SURFACE TEMPERATURE448 CALL get_var_phys( 'CDSW',qs) !--- SOIL MOISTURE364 CALL get_var_phys(tsrfvar,ts) !--- SURFACE TEMPERATURE 365 CALL get_var_phys(qsolvar,qs) !--- SOIL MOISTURE 449 366 CALL flinclo(fid_phys) 450 367 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) … … 463 380 ! 464 381 !------------------------------------------------------------------------------- 465 USE conf_dat_m, ONLY: conf_dat2d466 382 IMPLICIT NONE 467 383 !------------------------------------------------------------------------------- … … 474 390 !------------------------------------------------------------------------------- 475 391 SELECT CASE(title) 476 CASE( 'SP');tllm=0477 CASE( 'ST','CDSW'); tllm=llm_phys392 CASE(psrfvar); tllm=0 393 CASE(tsrfvar,qsolvar); tllm=llm_phys 478 394 END SELECT 479 395 IF(ALLOCATED(field)) RETURN 480 396 ALLOCATE(field(iml,jml)); field(:,:)=0. 481 397 CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana) 482 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)483 CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,&484 lon_in, lat_in, lon_in2, lat_in2,field)398 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 399 CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, & 400 lon_in, lat_in, field) 485 401 486 402 END SUBROUTINE get_var_phys … … 495 411 !------------------------------------------------------------------------------- 496 412 ! 497 SUBROUTINE interp_startvar(nam,ib ar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)413 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo) 498 414 ! 499 415 !------------------------------------------------------------------------------- 500 416 USE inter_barxy_m, ONLY: inter_barxy 501 USE grid_atob_m, ONLY: grille_m502 417 IMPLICIT NONE 503 418 !------------------------------------------------------------------------------- 504 419 ! Arguments: 505 420 CHARACTER(LEN=*), INTENT(IN) :: nam 506 LOGICAL, INTENT(IN) :: ib ar, ibeg421 LOGICAL, INTENT(IN) :: ibeg 507 422 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 508 423 REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) 509 REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1)510 424 REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) 511 425 REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) 512 426 !------------------------------------------------------------------------------- 513 427 ! Local variables: 514 CHARACTER(LEN=256) :: modname ="interp_startvar"428 CHARACTER(LEN=256) :: modname 515 429 INTEGER :: ii, jj, i1, j1, j2 516 430 REAL, ALLOCATABLE :: vtmp(:,:) 517 431 !------------------------------------------------------------------------------- 518 ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")519 jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")520 i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")521 j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1")522 j 2=SIZE(lat2)432 modname="interp_startvar" 433 ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") 434 jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") 435 i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") 436 j1=SIZE(varo,2); j2=SIZE(lat2) 523 437 ALLOCATE(vtmp(i1-1,j1)) 524 IF(ibar) THEN 525 IF(ibeg.AND.prt_level>1) THEN 526 WRITE(lunout,*)"--------------------------------------------------------" 527 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 528 WRITE(lunout,*)"--------------------------------------------------------" 529 END IF 530 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 531 ELSE 532 CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) 438 IF(ibeg.AND.prt_level>1) THEN 439 WRITE(lunout,*)"--------------------------------------------------------" 440 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 441 WRITE(lunout,*)"--------------------------------------------------------" 533 442 END IF 443 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 534 444 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 535 445
Note: See TracChangeset
for help on using the changeset viewer.