Changeset 2293 for LMDZ5/trunk/libf/phylmd/limit_netcdf.F90
- Timestamp:
- Jun 5, 2015, 9:16:07 PM (9 years ago)
- File:
-
- 1 edited
-
LMDZ5/trunk/libf/phylmd/limit_netcdf.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/limit_netcdf.F90
r2239 r2293 1 !2 ! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $3 !-------------------------------------------------------------------------------4 !5 1 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque) 6 #ifndef CPP_1D7 2 ! 8 3 !------------------------------------------------------------------------------- … … 21 16 ! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs) 22 17 !------------------------------------------------------------------------------- 18 #ifndef CPP_1D 23 19 USE control_mod 24 20 USE indice_sol_mod … … 32 28 NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT 33 29 USE inter_barxy_m, only: inter_barxy 34 use netcdf95, only: nf95_def_var, nf95_put_att, nf95_put_var 30 USE netcdf95, ONLY: nf95_def_var, nf95_put_att, nf95_put_var 31 USE grid_atob_m, ONLY: grille_m, rugosite, sea_ice 35 32 #endif 36 33 IMPLICIT NONE 37 34 !------------------------------------------------------------------------------- 38 35 ! Arguments: 39 #include "dimensions.h"40 #include "paramet.h"41 #include "iniprint.h"36 include "dimensions.h" 37 include "paramet.h" 38 include "iniprint.h" 42 39 LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation 43 40 LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag … … 49 46 !------------------------------------------------------------------------------- 50 47 ! Local variables: 51 #include "logic.h"52 #include "comvert.h"53 #include "comgeom2.h"54 #include "comconst.h"48 include "logic.h" 49 include "comvert.h" 50 include "comgeom2.h" 51 include "comconst.h" 55 52 56 53 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 57 CHARACTER(LEN=25) :: icefile, sstfile, dumstr 58 CHARACTER(LEN=25), PARAMETER :: famipsst ='amipbc_sst_1x1.nc ', & 59 famipsic ='amipbc_sic_1x1.nc ', & 60 fcpldsst ='cpl_atm_sst.nc ', & 61 fcpldsic ='cpl_atm_sic.nc ', & 62 fhistsst ='histmth_sst.nc ', & 63 fhistsic ='histmth_sic.nc ', & 64 frugo ='Rugos.nc ', & 65 falbe ='Albedo.nc ', & 66 feraisstk='sstk.nc ', & 67 feraici ='ci.nc ' 54 CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam 55 CHARACTER(LEN=20), PARAMETER :: & 56 fsst(4)=['amipbc_sst_1x1.nc ','cpl_atm_sst.nc ','histmth_sst.nc '& 57 ,'sstk.nc '] 58 CHARACTER(LEN=20), PARAMETER :: & 59 fsic(4)=['amipbc_sic_1x1.nc ','cpl_atm_sic.nc ','histmth_sic.nc '& 60 ,'ci.nc '] 61 CHARACTER(LEN=10), PARAMETER :: & 62 vsst(4)=['tosbcs ','SISUTESW ','tsol_oce ','sstk '], & 63 vsic(4)=['sicbcs ','SIICECOV ','pourc_sic ','ci '] 64 CHARACTER(LEN=20), PARAMETER :: frugo='Rugos.nc ', & 65 falbe='Albedo.nc ' 66 CHARACTER(LEN=10), PARAMETER :: vrug='RUGOS ', valb='ALBEDO ' 68 67 CHARACTER(LEN=10) :: varname 69 68 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 70 REAL, DIMENSION(klon) :: fi_ice, verif 71 REAL, DIMENSION(:,:), POINTER :: phy_rug=>NULL(), phy_ice=>NULL() 72 REAL, DIMENSION(:,:), POINTER :: phy_sst=>NULL(), phy_alb=>NULL() 73 REAL, DIMENSION(:,:), ALLOCATABLE :: phy_bil 74 REAL, DIMENSION(:,:,:), ALLOCATABLE :: pctsrf_t 75 INTEGER :: nbad 69 REAL :: fi_ice(klon), verif(klon) 70 REAL, POINTER :: phy_rug(:,:)=>NULL(), phy_ice(:,:)=>NULL() 71 REAL, POINTER :: phy_sst(:,:)=>NULL(), phy_alb(:,:)=>NULL() 72 REAL, ALLOCATABLE :: phy_bil(:,:), pctsrf_t(:,:,:) 73 INTEGER :: nbad 76 74 77 75 !--- VARIABLES FOR OUTPUT FILE WRITING ----------------------------------------- 78 INTEGER :: ierr, nid, ndim, ntim, k 79 INTEGER, DIMENSION(2) :: dims 76 INTEGER :: ierr, nid, ndim, ntim, k, dims(2), ix_sic, ix_sst 80 77 INTEGER :: id_tim, id_SST, id_BILS, id_RUG, id_ALB 81 78 INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC, varid_longitude, varid_latitude … … 89 86 NF90_FORMAT=NF90_FLOAT 90 87 #endif 91 92 pi = 4.*ATAN(1.)93 rad = 6371229.94 daysec= 86400.95 omeg = 2.*pi/daysec96 g = 9.897 kappa = 0.285714398 cpp = 1004.7088599 dtvr = daysec/REAL(day_step)100 88 CALL inigeom 101 89 … … 104 92 105 93 !--- RUGOSITY TREATMENT -------------------------------------------------------- 106 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' 107 varname='RUGOS' 108 CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 94 CALL msg(1,'Traitement de la rugosite') 95 CALL get_2Dfield(frugo,vrug,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 109 96 110 97 !--- OCEAN TREATMENT ----------------------------------------------------------- 111 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'98 CALL msg(1,'Traitement de la glace oceanique') 112 99 113 100 ! Input SIC file selection 114 101 ! Open file only to test if available 115 IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 116 icefile=TRIM(famipsic) 117 varname='sicbcs' 118 ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 119 icefile=TRIM(fcpldsic) 120 varname='SIICECOV' 121 ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 122 icefile=TRIM(fhistsic) 123 varname='pourc_sic' 124 ELSE IF ( NF90_OPEN(TRIM(feraici),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 125 icefile=TRIM(feraici) 126 varname='ci' 127 ELSE 102 DO ix_sic=1,SIZE(fsic) 103 IF ( NF90_OPEN(TRIM(fsic(ix_sic)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 104 icefile=fsic(ix_sic); varname=vsic(ix_sic); EXIT 105 END IF 106 END DO 107 IF(ix_sic==SIZE(fsic)+1) THEN 128 108 WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' 129 WRITE(lunout,*) 'One of following files must be avail ible : ',trim(famipsic),', ',trim(fcpldsic),', ', &130 trim(fhistsic), trim(feraici)109 WRITE(lunout,*) 'One of following files must be available : ' 110 DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO 131 111 CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) 132 112 END IF 133 ierr=NF90_CLOSE(nid)134 IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)113 CALL ncerr(NF90_CLOSE(nid),icefile) 114 CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile)) 135 115 136 116 CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) … … 143 123 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 144 124 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 145 IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER) 125 SELECT CASE(ix_sic) 126 CASE(2) ! SIC=pICE*(1-LIC-TER) (CPL) 146 127 pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 147 ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE128 CASE(3) ! SIC=pICE (HIST) 148 129 pctsrf_t(:,is_sic,k)=fi_ice(:) 149 ELSE ! icefile==famipsic ! SIC=pICE-LIC130 CASE DEFAULT ! SIC=pICE-LIC (AMIP,ERAI) 150 131 pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k) 151 END IF132 END SELECT 152 133 WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0. 153 134 WHERE(1.0-zmasq<EPSFRA) … … 167 148 END WHERE 168 149 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 169 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad170 nbad=COUNT( abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)150 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb points = ',nbad 151 nbad=COUNT(ABS(SUM(pctsrf_t(:,:,k),DIM=2)-1.0)>EPSFRA) 171 152 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 172 153 END DO … … 174 155 175 156 !--- SST TREATMENT ------------------------------------------------------------- 176 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'157 CALL msg(1,'Traitement de la sst') 177 158 178 159 ! Input SST file selection 179 160 ! Open file only to test if available 180 IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 181 sstfile=TRIM(famipsst) 182 varname='tosbcs' 183 ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 184 sstfile=TRIM(fcpldsst) 185 varname='SISUTESW' 186 ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 187 sstfile=TRIM(fhistsst) 188 varname='tsol_oce' 189 ELSE IF ( NF90_OPEN(TRIM(feraisstk),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 190 sstfile=TRIM(feraisstk) 191 varname='sstk' 192 ELSE 161 DO ix_sst=1,SIZE(fsst) 162 IF ( NF90_OPEN(TRIM(fsst(ix_sst)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 163 sstfile=fsst(ix_sst); varname=vsst(ix_sst); EXIT 164 END IF 165 END DO 166 IF(ix_sst==SIZE(fsst)+1) THEN 193 167 WRITE(lunout,*) 'ERROR! No sst input file was found.' 194 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst),trim(feraisstk) 168 WRITE(lunout,*) 'One of following files must be available : ' 169 DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO 195 170 CALL abort_gcm('limit_netcdf','No sst file was found',1) 196 171 END IF 197 ierr=NF90_CLOSE(nid)198 IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)172 CALL ncerr(NF90_CLOSE(nid),sstfile) 173 CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile)) 199 174 200 175 CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) 201 176 202 177 !--- ALBEDO TREATMENT ---------------------------------------------------------- 203 IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' 204 varname='ALBEDO' 205 CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) 178 CALL msg(1,'Traitement de l albedo') 179 CALL get_2Dfield(falbe,valb,'ALB',interbar,ndays,phy_alb) 206 180 207 181 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- … … 209 183 210 184 !--- OUTPUT FILE WRITING ------------------------------------------------------- 211 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' 185 CALL msg(5,'Ecriture du fichier limit : debut') 186 fnam="limit.nc" 212 187 213 188 !--- File creation 214 ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)215 ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")189 CALL ncerr(NF90_CREATE(fnam,NF90_CLOBBER,nid),fnam) 190 CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam) 216 191 217 192 !--- Dimensions creation 218 ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)219 ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)220 221 dims= (/ndim,ntim/)193 CALL ncerr(NF90_DEF_DIM(nid,"points_physiques",klon,ndim),fnam) 194 CALL ncerr(NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim),fnam) 195 196 dims=[ndim,ntim] 222 197 223 198 !--- Variables creation 224 ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,(/ntim/),id_tim)225 ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)226 ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)227 ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)228 ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)229 ierr=NF90_DEF_VAR(nid,"SST", NF90_FORMAT,dims,id_SST)230 ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)231 ierr=NF90_DEF_VAR(nid,"ALB", NF90_FORMAT,dims,id_ALB)232 ierr=NF90_DEF_VAR(nid,"RUG", NF90_FORMAT,dims,id_RUG)199 CALL ncerr(NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,[ntim],id_tim),fnam) 200 CALL ncerr(NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE),fnam) 201 CALL ncerr(NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC),fnam) 202 CALL ncerr(NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER),fnam) 203 CALL ncerr(NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC),fnam) 204 CALL ncerr(NF90_DEF_VAR(nid,"SST", NF90_FORMAT,dims,id_SST),fnam) 205 CALL ncerr(NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS),fnam) 206 CALL ncerr(NF90_DEF_VAR(nid,"ALB", NF90_FORMAT,dims,id_ALB),fnam) 207 CALL ncerr(NF90_DEF_VAR(nid,"RUG", NF90_FORMAT,dims,id_RUG),fnam) 233 208 call nf95_def_var(nid, "longitude", NF90_FLOAT, ndim, varid_longitude) 234 call nf95_def_var(nid, "latitude", NF90_FLOAT, ndim, varid_latitude)209 call nf95_def_var(nid, "latitude", NF90_FLOAT, ndim, varid_latitude) 235 210 236 211 !--- Attributes creation 237 ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")238 ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")239 ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")240 ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")241 ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")242 ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")243 ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")244 ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")245 ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")212 CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam) 213 CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam) 214 CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam) 215 CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam) 216 CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam) 217 CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam) 218 CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam) 219 CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam) 220 CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam) 246 221 247 222 call nf95_put_att(nid, varid_longitude, "standard_name", "longitude") … … 251 226 call nf95_put_att(nid, varid_latitude, "units", "degrees_north") 252 227 253 ierr=NF90_ENDDEF(nid)228 CALL ncerr(NF90_ENDDEF(nid),fnam) 254 229 255 230 !--- Variables saving 256 ierr=NF90_PUT_VAR(nid,id_tim,(/(REAL(k),k=1,ndays)/))257 ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))258 ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))259 ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))260 ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))261 ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))262 ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))263 ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))264 ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))231 CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam) 232 CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam) 233 CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam) 234 CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam) 235 CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam) 236 CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam) 237 CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam) 238 CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam) 239 CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam) 265 240 call nf95_put_var(nid, varid_longitude, rlon) 266 241 call nf95_put_var(nid, varid_latitude, rlat) 267 242 268 ierr=NF90_CLOSE(nid)269 270 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'243 CALL ncerr(NF90_CLOSE(nid),fnam) 244 245 CALL msg(6,'Ecriture du fichier limit : fin') 271 246 272 247 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) … … 296 271 USE dimphy, ONLY : klon 297 272 USE phys_state_var_mod, ONLY : pctsrf 273 USE conf_dat_m, ONLY: conf_dat2d 298 274 USE control_mod 299 275 USE pchsp_95_m, only: pchsp_95 … … 303 279 304 280 IMPLICIT NONE 305 #include "dimensions.h"306 #include "paramet.h"307 #include "comgeom2.h"308 #include "iniprint.h"281 include "dimensions.h" 282 include "paramet.h" 283 include "comgeom2.h" 284 include "iniprint.h" 309 285 !----------------------------------------------------------------------------- 310 286 ! Arguments: … … 320 296 ! Local variables: 321 297 !--- NetCDF 322 INTEGER :: ncid, varid! NetCDF identifiers323 CHARACTER(LEN=30) :: dnam! dimension name298 INTEGER :: ncid, varid ! NetCDF identifiers 299 CHARACTER(LEN=30) :: dnam ! dimension name 324 300 !--- dimensions 325 INTEGER , DIMENSION(4) :: dids! NetCDF dimensions identifiers326 REAL, ALLOCATABLE , DIMENSION(:) :: dlon_ini! initial longitudes vector327 REAL, ALLOCATABLE , DIMENSION(:) :: dlat_ini! initial latitudes vector328 REAL, POINTER , DIMENSION(:) :: dlon, dlat! reordered lon/lat vectors301 INTEGER :: dids(4) ! NetCDF dimensions identifiers 302 REAL, ALLOCATABLE :: dlon_ini(:) ! initial longitudes vector 303 REAL, ALLOCATABLE :: dlat_ini(:) ! initial latitudes vector 304 REAL, POINTER :: dlon(:), dlat(:) ! reordered lon/lat vectors 329 305 !--- fields 330 INTEGER :: imdep, jmdep, lmdep 331 REAL, ALLOCATABLE , DIMENSION(:, :) :: champ! wanted field on initial grid332 REAL, ALLOCATABLE , DIMENSION(:) :: yder, timeyear333 REAL , DIMENSION(iim, jjp1) :: champint! interpolated field334 REAL, ALLOCATABLE , DIMENSION(:, :, :) :: champtime335 REAL, ALLOCATABLE , DIMENSION(:, :, :) :: champan306 INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ' 307 REAL, ALLOCATABLE :: champ(:,:) ! wanted field on initial grid 308 REAL, ALLOCATABLE :: yder(:), timeyear(:) 309 REAL :: champint(iim,jjp1) ! interpolated field 310 REAL, ALLOCATABLE :: champtime(:,:,:) 311 REAL, ALLOCATABLE :: champan(:,:,:) 336 312 !--- input files 337 CHARACTER(LEN=20) :: cal_in! calendar338 CHARACTER(LEN=20) :: unit_sic! attribute unit in sea-ice file339 INTEGER :: ndays_in! number of days313 CHARACTER(LEN=20) :: cal_in ! calendar 314 CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file 315 INTEGER :: ndays_in ! number of days 340 316 !--- misc 341 INTEGER :: i, j, k, l! loop counters342 REAL, ALLOCATABLE , DIMENSION(:, :) :: work! used for extrapolation343 CHARACTER(LEN=25) :: title! for messages344 LOGICAL :: extrp! flag for extrapolation345 LOGICAL :: oldice! flag for old way ice computation346 REAL 317 INTEGER :: i, j, k, l ! loop counters 318 REAL, ALLOCATABLE :: work(:,:) ! used for extrapolation 319 CHARACTER(LEN=25) :: title ! for messages 320 LOGICAL :: extrp ! flag for extrapolation 321 LOGICAL :: oldice ! flag for old way ice computation 322 REAL :: chmin, chmax 347 323 INTEGER ierr 348 324 integer n_extrap ! number of extrapolated points … … 369 345 370 346 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- 371 IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam372 ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr,fnam)373 ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr,fnam)374 ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr,fnam)347 CALL msg(5,' Now reading file : '//TRIM(fnam)) 348 CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam) 349 CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam) 350 CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam) 375 351 376 352 !--- Read unit for sea-ice variable only 377 353 IF (mode=='SIC') THEN 378 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) 379 IF(ierr/=NF90_NOERR) THEN 380 IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' 381 unit_sic='X' 382 ELSE 383 IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic 384 END IF 354 IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN 355 CALL msg(5,'No unit in sea-ice file. Take percentage as default value') 356 unit_sic='X' 357 ELSE 358 CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic)) 359 END IF 385 360 END IF 386 361 387 362 !--- Longitude 388 ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)389 CALL ncerr(ierr, fnam);ALLOCATE(dlon_ini(imdep), dlon(imdep))390 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)391 ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam)392 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep363 CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam) 364 ALLOCATE(dlon_ini(imdep), dlon(imdep)) 365 CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam) 366 CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam) 367 CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep) 393 368 394 369 !--- Latitude 395 ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)396 CALL ncerr(ierr, fnam);ALLOCATE(dlat_ini(jmdep), dlat(jmdep))397 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)398 ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam)399 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep370 CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam) 371 ALLOCATE(dlat_ini(jmdep), dlat(jmdep)) 372 CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam) 373 CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam) 374 CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep) 400 375 401 376 !--- Time (variable is not needed - it is rebuilt - but calendar is) 402 ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)403 CALL ncerr(ierr, fnam);ALLOCATE(timeyear(lmdep))404 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)377 CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam) 378 ALLOCATE(timeyear(lmdep)) 379 CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam) 405 380 cal_in=' ' 406 ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in) 407 IF(ierr/=NF90_NOERR) THEN 381 IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN 408 382 SELECT CASE(mode) 409 383 CASE('RUG', 'ALB'); cal_in='360d' 410 384 CASE('SIC', 'SST'); cal_in='gregorian' 411 385 END SELECT 412 IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & 413 // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' 414 END IF 415 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & 416 cal_in 417 386 CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '& 387 &//TRIM(fnam)//'. Choosing default value.') 388 END IF 389 CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep) 418 390 419 391 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- … … 430 402 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep)) 431 403 IF(extrp) ALLOCATE(work(imdep, jmdep)) 432 433 IF (prt_level>5) WRITE(lunout, *) 434 IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' 435 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) 404 CALL msg(5,'') 405 CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.') 406 CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam) 436 407 DO l=1, lmdep 437 ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/)) 438 CALL ncerr(ierr, fnam) 439 CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, & 440 champ, ibar) 441 442 IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, & 443 work) 408 CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam) 409 CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar) 410 IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work) 444 411 445 412 IF(ibar .AND. .NOT.oldice) THEN 446 IF(l==1 .AND. prt_level>5) THEN447 WRITE(lunout, *) '-------------------------------------------------------------------------'448 WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'449 WRITE(lunout, *) '-------------------------------------------------------------------------'413 IF(l==1) THEN 414 CALL msg(5,"----------------------------------------------------------") 415 CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$") 416 CALL msg(5,"----------------------------------------------------------") 450 417 END IF 451 418 IF(mode=='RUG') champ=LOG(champ) 452 CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), & 453 rlatv, champint) 419 CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint) 454 420 IF(mode=='RUG') THEN 455 421 champint=EXP(champint) … … 458 424 ELSE 459 425 SELECT CASE(mode) 460 CASE('RUG'); CALL rugosite(imdep, jmdep, dlon, dlat, champ, & 461 iim, jjp1, rlonv, rlatu, champint, mask) 462 CASE('SIC'); CALL sea_ice (imdep, jmdep, dlon, dlat, champ, & 463 iim, jjp1, rlonv, rlatu, champint) 464 CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, & 465 iim, jjp1, rlonv, rlatu, champint) 426 CASE('RUG'); CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask) 427 CASE('SIC'); CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint) 428 CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint) 466 429 END SELECT 467 430 END IF 468 431 champtime(:, :, l)=champint 469 432 END DO 470 ierr=NF90_CLOSE(ncid); CALL ncerr(ierr, fnam)433 CALL ncerr(NF90_CLOSE(ncid), fnam) 471 434 472 435 DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ) … … 474 437 475 438 !--- TIME INTERPOLATION ------------------------------------------------------ 476 IF 439 IF(prt_level>5) THEN 477 440 WRITE(lunout, *) 478 441 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' … … 508 471 !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- 509 472 IF(mode=='SST') THEN 510 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'473 CALL msg(5,'Filtrage de la SST: SST >= 271.38') 511 474 WHERE(champan<271.38) champan=271.38 512 475 END IF … … 514 477 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ------------------------------------- 515 478 IF(mode=='SIC') THEN 516 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'479 CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0') 517 480 518 481 IF (unit_sic=='1') THEN 519 482 ! Nothing to be done for sea-ice field is already in fraction of 1 520 483 ! This is the case for sea-ice in file cpl_atm_sic.nc 521 IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'484 CALL msg(5,'Sea-ice field already in fraction of 1') 522 485 ELSE 523 486 ! Convert sea ice from percentage to fraction of 1 524 IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.'487 CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.') 525 488 champan(:, :, :)=champan(:, :, :)/100. 526 489 END IF … … 541 504 ! 542 505 !------------------------------------------------------------------------------- 543 544 506 545 507 … … 620 582 621 583 ELSE 622 mnth= (/(m,m=1,nm,nd/nm)/)584 mnth=[(m,m=1,nm,nd/nm)] 623 585 END IF 624 586 … … 634 596 635 597 598 599 !------------------------------------------------------------------------------- 600 ! 601 SUBROUTINE msg(lev,str1,i,str2) 602 ! 603 !------------------------------------------------------------------------------- 604 ! Arguments: 605 INTEGER, INTENT(IN) :: lev 606 CHARACTER(LEN=*), INTENT(IN) :: str1 607 INTEGER, OPTIONAL, INTENT(IN) :: i 608 CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2 609 !------------------------------------------------------------------------------- 610 IF(prt_level>lev) THEN 611 IF(PRESENT(str2)) THEN 612 WRITE(lunout,*) TRIM(str1), i, TRIM(str2) 613 ELSE IF(PRESENT(i)) THEN 614 WRITE(lunout,*) TRIM(str1), i 615 ELSE 616 WRITE(lunout,*) TRIM(str1) 617 END IF 618 END IF 619 620 END SUBROUTINE msg 621 ! 622 !------------------------------------------------------------------------------- 623 624 636 625 !------------------------------------------------------------------------------- 637 626 !
Note: See TracChangeset
for help on using the changeset viewer.
