[38] | 1 | subroutine writediagsoil(ngrid,name,title,units,dimpx,px) |
---|
| 2 | |
---|
| 3 | ! Write variable 'name' to NetCDF file 'diagsoil.nc'. |
---|
| 4 | ! The variable may be 3D (lon,lat,depth) subterranean field, |
---|
| 5 | ! a 2D (lon,lat) surface field, or a simple scalar (0D variable). |
---|
| 6 | ! |
---|
| 7 | ! Calls to 'writediagsoil' can originate from anywhere in the program; |
---|
| 8 | ! An initialisation of variable 'name' is done if it is the first time |
---|
| 9 | ! that this routine is called with given 'name'; otherwise data is appended |
---|
| 10 | ! (yielding the sought time series of the variable) |
---|
| 11 | |
---|
| 12 | ! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4 |
---|
| 13 | |
---|
[1047] | 14 | use comsoil_h, only: nsoilmx |
---|
[1130] | 15 | use control_mod, only: ecritphy, day_step, iphysiq |
---|
| 16 | use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather |
---|
| 17 | use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo |
---|
[1047] | 18 | |
---|
[38] | 19 | implicit none |
---|
| 20 | |
---|
| 21 | #include"dimensions.h" |
---|
| 22 | #include"paramet.h" |
---|
[1130] | 23 | !#include"control.h" |
---|
[1047] | 24 | !#include"comsoil.h" |
---|
[38] | 25 | #include"netcdf.inc" |
---|
| 26 | |
---|
| 27 | ! Arguments: |
---|
| 28 | integer,intent(in) :: ngrid ! number of (horizontal) points of physics grid |
---|
| 29 | ! i.e. ngrid = 2+(jjm-1)*iim - 1/jjm |
---|
| 30 | character(len=*),intent(in) :: name ! 'name' of the variable |
---|
| 31 | character(len=*),intent(in) :: title ! 'long_name' attribute of the variable |
---|
| 32 | character(len=*),intent(in) :: units ! 'units' attribute of the variable |
---|
| 33 | integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0) |
---|
| 34 | real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable |
---|
[1047] | 35 | ! Note: nsoilmx is a parameter set in 'comsoil_h' |
---|
[38] | 36 | |
---|
| 37 | ! Local variables: |
---|
| 38 | real*4,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data |
---|
[1047] | 39 | ! Note iip1,jjp1 known from paramet.h; nsoilmx known from comsoil_h |
---|
[38] | 40 | real*4,dimension(iip1,jjp1) :: data2 ! to store 2D data |
---|
| 41 | real*4 :: data0 ! to store 0D data |
---|
| 42 | integer :: i,j,l ! for loops |
---|
| 43 | integer :: ig0 |
---|
| 44 | |
---|
| 45 | real*4,save :: date ! time counter (in elapsed days) |
---|
| 46 | integer,save :: isample ! sample rate at which data is to be written to output |
---|
| 47 | integer,save :: ntime=0 ! counter to internally store time steps |
---|
| 48 | character(len=20),save :: firstname="1234567890" |
---|
| 49 | integer,save :: zitau=0 |
---|
| 50 | |
---|
| 51 | character(len=30) :: filename="diagsoil.nc" |
---|
| 52 | |
---|
| 53 | ! NetCDF stuff: |
---|
| 54 | integer :: nid ! NetCDF output file ID |
---|
| 55 | integer :: varid ! NetCDF ID of a variable |
---|
| 56 | integer :: ierr ! NetCDF routines return code |
---|
| 57 | integer,dimension(4) :: id ! NetCDF IDs of the dimensions of the variable |
---|
| 58 | integer,dimension(4) :: edges,corners |
---|
| 59 | |
---|
[1130] | 60 | #ifdef CPP_PARA |
---|
| 61 | ! Added to work in parallel mode |
---|
| 62 | real dx3_glop(klon_glo,nsoilmx) |
---|
| 63 | real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set |
---|
| 64 | real dx2_glop(klon_glo) |
---|
| 65 | real dx2_glo(iim,jjp1) ! to store a global 2D (surface) data set |
---|
| 66 | real px2(ngrid) |
---|
| 67 | #endif |
---|
| 68 | |
---|
[38] | 69 | ! 1. Initialization step |
---|
| 70 | if (firstname.eq."1234567890") then |
---|
| 71 | ! Store 'name' as 'firstname' |
---|
| 72 | firstname=name |
---|
| 73 | ! From now on, if 'name'.eq.'firstname', then it is a new time cycle |
---|
| 74 | |
---|
| 75 | ! just to be sure, check that firstnom is large enough to hold nom |
---|
| 76 | if (len_trim(firstname).lt.len_trim(name)) then |
---|
| 77 | write(*,*) "writediagsoil: Error !!!" |
---|
| 78 | write(*,*) " firstname string not long enough!!" |
---|
| 79 | write(*,*) " increase its size to at least ",len_trim(name) |
---|
| 80 | stop |
---|
| 81 | endif |
---|
| 82 | |
---|
| 83 | ! Set output sample rate |
---|
| 84 | isample=int(ecritphy) ! same as for diagfi outputs |
---|
| 85 | ! Note ecritphy is known from control.h |
---|
| 86 | |
---|
| 87 | ! Create output NetCDF file |
---|
[1130] | 88 | if (is_master) then |
---|
| 89 | ierr=NF_CREATE(filename,NF_CLOBBER,nid) |
---|
| 90 | if (ierr.ne.NF_NOERR) then |
---|
[38] | 91 | write(*,*)'writediagsoil: Error, failed creating file '//trim(filename) |
---|
| 92 | stop |
---|
[1130] | 93 | endif |
---|
| 94 | endif ! of if (is_master) |
---|
| 95 | |
---|
[38] | 96 | ! Define dimensions and axis attributes |
---|
[1047] | 97 | call iniwritesoil(nid,ngrid) |
---|
[38] | 98 | |
---|
| 99 | ! set zitau to -1 to be compatible with zitau incrementation step below |
---|
| 100 | zitau=-1 |
---|
| 101 | |
---|
| 102 | else |
---|
| 103 | ! If not an initialization call, simply open the NetCDF file |
---|
[1130] | 104 | if (is_master) then |
---|
| 105 | ierr=NF_OPEN(filename,NF_WRITE,nid) |
---|
| 106 | endif |
---|
[38] | 107 | endif ! of if (firstname.eq."1234567890") |
---|
| 108 | |
---|
| 109 | ! 2. Increment local time counter, if necessary |
---|
| 110 | if (name.eq.firstname) then |
---|
| 111 | ! if we run across 'firstname', then it is a new time step |
---|
| 112 | zitau=zitau+iphysiq |
---|
| 113 | ! Note iphysiq is known from control.h |
---|
| 114 | endif |
---|
| 115 | |
---|
| 116 | ! 3. Write data, if the time index matches the sample rate |
---|
| 117 | if (mod(zitau+1,isample).eq.0) then |
---|
| 118 | |
---|
| 119 | ! 3.1 If first call at this date, update 'time' variable |
---|
| 120 | if (name.eq.firstname) then |
---|
| 121 | ntime=ntime+1 |
---|
[1130] | 122 | date=float(zitau+1)/float(day_step) |
---|
[38] | 123 | ! Note: day_step is known from control.h |
---|
| 124 | |
---|
[1130] | 125 | if (is_master) then |
---|
| 126 | ! Get NetCDF ID for "time" |
---|
| 127 | ierr=NF_INQ_VARID(nid,"time",varid) |
---|
| 128 | ! Add the current value of date to the "time" array |
---|
[38] | 129 | !#ifdef NC_DOUBLE |
---|
[1130] | 130 | ! ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date) |
---|
[38] | 131 | !#else |
---|
[1130] | 132 | ierr=NF_PUT_VARA_REAL(nid,varid,ntime,1,date) |
---|
[38] | 133 | !#endif |
---|
[1130] | 134 | if (ierr.ne.NF_NOERR) then |
---|
[38] | 135 | write(*,*)"writediagsoil: Failed writing date to time variable" |
---|
| 136 | stop |
---|
[1130] | 137 | endif |
---|
| 138 | endif ! of if (is_master) |
---|
[38] | 139 | endif ! of if (name.eq.firstname) |
---|
| 140 | |
---|
| 141 | ! 3.2 Write the variable to the NetCDF file |
---|
| 142 | if (dimpx.eq.3) then ! Case of a 3D variable |
---|
| 143 | ! A. Recast data along 'dynamics' grid |
---|
[1130] | 144 | #ifdef CPP_PARA |
---|
| 145 | ! gather field on a "global" (without redundant longitude) array |
---|
| 146 | call Gather(px,dx3_glop) |
---|
| 147 | !$OMP MASTER |
---|
| 148 | if (is_mpi_root) then |
---|
| 149 | call Grid1Dto2D_glo(dx3_glop,dx3_glo) |
---|
| 150 | ! copy dx3_glo() to dx3(:) and add redundant longitude |
---|
| 151 | data3(1:iim,:,:)=dx3_glo(1:iim,:,:) |
---|
| 152 | data3(iip1,:,:)=data3(1,:,:) |
---|
| 153 | endif |
---|
| 154 | !$OMP END MASTER |
---|
| 155 | !$OMP BARRIER |
---|
| 156 | #else |
---|
[38] | 157 | do l=1,nsoilmx |
---|
| 158 | ! handle the poles |
---|
| 159 | do i=1,iip1 |
---|
| 160 | data3(i,1,l)=px(1,l) |
---|
| 161 | data3(i,jjp1,l)=px(ngrid,l) |
---|
| 162 | enddo |
---|
| 163 | ! rest of the grid |
---|
| 164 | do j=2,jjm |
---|
| 165 | ig0=1+(j-2)*iim |
---|
| 166 | do i=1,iim |
---|
| 167 | data3(i,j,l)=px(ig0+i,l) |
---|
| 168 | enddo |
---|
| 169 | data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude |
---|
| 170 | enddo |
---|
| 171 | enddo |
---|
[1130] | 172 | #endif |
---|
[38] | 173 | |
---|
| 174 | ! B. Write (append) the variable to the NetCDF file |
---|
[1130] | 175 | if (is_master) then |
---|
[38] | 176 | ! B.1. Get the ID of the variable |
---|
| 177 | ierr=NF_INQ_VARID(nid,name,varid) |
---|
| 178 | if (ierr.ne.NF_NOERR) then |
---|
| 179 | ! If we failed geting the variable's ID, we assume it is because |
---|
| 180 | ! the variable doesn't exist yet and must be created. |
---|
| 181 | ! Start by obtaining corresponding dimensions IDs |
---|
| 182 | ierr=NF_INQ_DIMID(nid,"longitude",id(1)) |
---|
| 183 | ierr=NF_INQ_DIMID(nid,"latitude",id(2)) |
---|
| 184 | ierr=NF_INQ_DIMID(nid,"depth",id(3)) |
---|
| 185 | ierr=NF_INQ_DIMID(nid,"time",id(4)) |
---|
| 186 | ! Tell the world about it |
---|
| 187 | write(*,*) "=====================" |
---|
| 188 | write(*,*) "writediagsoil: creating variable "//trim(name) |
---|
| 189 | call def_var(nid,name,title,units,4,id,varid,ierr) |
---|
| 190 | endif ! of if (ierr.ne.NF_NOERR) |
---|
| 191 | |
---|
| 192 | ! B.2. Prepare things to be able to write/append the variable |
---|
| 193 | corners(1)=1 |
---|
| 194 | corners(2)=1 |
---|
| 195 | corners(3)=1 |
---|
| 196 | corners(4)=ntime |
---|
| 197 | |
---|
| 198 | edges(1)=iip1 |
---|
| 199 | edges(2)=jjp1 |
---|
| 200 | edges(3)=nsoilmx |
---|
| 201 | edges(4)=1 |
---|
| 202 | |
---|
| 203 | ! B.3. Write the slab of data |
---|
| 204 | !#ifdef NC_DOUBLE |
---|
| 205 | ! ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3) |
---|
| 206 | !#else |
---|
| 207 | ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3) |
---|
| 208 | !#endif |
---|
| 209 | if (ierr.ne.NF_NOERR) then |
---|
| 210 | write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//& |
---|
| 211 | " to file "//trim(filename)//" at time",date |
---|
| 212 | endif |
---|
[1130] | 213 | endif ! of if (is_master) |
---|
[38] | 214 | |
---|
| 215 | elseif (dimpx.eq.2) then ! Case of a 2D variable |
---|
[1130] | 216 | |
---|
[38] | 217 | ! A. Recast data along 'dynamics' grid |
---|
[1130] | 218 | #ifdef CPP_PARA |
---|
| 219 | ! gather field on a "global" (without redundant longitude) array |
---|
| 220 | px2(:)=px(:,1) |
---|
| 221 | call Gather(px2,dx2_glop) |
---|
| 222 | !$OMP MASTER |
---|
| 223 | if (is_mpi_root) then |
---|
| 224 | call Grid1Dto2D_glo(dx2_glop,dx2_glo) |
---|
| 225 | ! copy dx3_glo() to dx3(:) and add redundant longitude |
---|
| 226 | data2(1:iim,:)=dx2_glo(1:iim,:) |
---|
| 227 | data2(iip1,:)=data2(1,:) |
---|
| 228 | endif |
---|
| 229 | !$OMP END MASTER |
---|
| 230 | !$OMP BARRIER |
---|
| 231 | #else |
---|
[38] | 232 | ! handle the poles |
---|
| 233 | do i=1,iip1 |
---|
| 234 | data2(i,1)=px(1,1) |
---|
| 235 | data2(i,jjp1)=px(ngrid,1) |
---|
| 236 | enddo |
---|
| 237 | ! rest of the grid |
---|
| 238 | do j=2,jjm |
---|
| 239 | ig0=1+(j-2)*iim |
---|
| 240 | do i=1,iim |
---|
| 241 | data2(i,j)=px(ig0+i,1) |
---|
| 242 | enddo |
---|
| 243 | data2(iip1,j)=data2(1,j) ! extra (modulo) longitude |
---|
| 244 | enddo |
---|
[1130] | 245 | #endif |
---|
[38] | 246 | |
---|
| 247 | ! B. Write (append) the variable to the NetCDF file |
---|
[1130] | 248 | if (is_master) then |
---|
[38] | 249 | ! B.1. Get the ID of the variable |
---|
| 250 | ierr=NF_INQ_VARID(nid,name,varid) |
---|
| 251 | if (ierr.ne.NF_NOERR) then |
---|
| 252 | ! If we failed geting the variable's ID, we assume it is because |
---|
| 253 | ! the variable doesn't exist yet and must be created. |
---|
| 254 | ! Start by obtaining corresponding dimensions IDs |
---|
| 255 | ierr=NF_INQ_DIMID(nid,"longitude",id(1)) |
---|
| 256 | ierr=NF_INQ_DIMID(nid,"latitude",id(2)) |
---|
| 257 | ierr=NF_INQ_DIMID(nid,"time",id(3)) |
---|
| 258 | ! Tell the world about it |
---|
| 259 | write(*,*) "=====================" |
---|
| 260 | write(*,*) "writediagsoil: creating variable "//trim(name) |
---|
| 261 | call def_var(nid,name,title,units,3,id,varid,ierr) |
---|
| 262 | endif ! of if (ierr.ne.NF_NOERR) |
---|
| 263 | |
---|
| 264 | ! B.2. Prepare things to be able to write/append the variable |
---|
| 265 | corners(1)=1 |
---|
| 266 | corners(2)=1 |
---|
| 267 | corners(3)=ntime |
---|
| 268 | |
---|
| 269 | edges(1)=iip1 |
---|
| 270 | edges(2)=jjp1 |
---|
| 271 | edges(3)=1 |
---|
| 272 | |
---|
| 273 | ! B.3. Write the slab of data |
---|
| 274 | !#ifdef NC_DOUBLE |
---|
| 275 | ! ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2) |
---|
| 276 | !#else |
---|
| 277 | ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2) |
---|
| 278 | !#endif |
---|
| 279 | if (ierr.ne.NF_NOERR) then |
---|
| 280 | write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//& |
---|
| 281 | " to file "//trim(filename)//" at time",date |
---|
| 282 | endif |
---|
[1130] | 283 | endif ! of if (is_master) |
---|
[38] | 284 | |
---|
| 285 | elseif (dimpx.eq.0) then ! Case of a 0D variable |
---|
[1130] | 286 | #ifdef CPP_PARA |
---|
| 287 | write(*,*) "writediagsoil: dimps==0 case not implemented in // mode!!" |
---|
| 288 | stop |
---|
| 289 | #endif |
---|
[38] | 290 | ! A. Copy data value |
---|
| 291 | data0=px(1,1) |
---|
| 292 | |
---|
| 293 | ! B. Write (append) the variable to the NetCDF file |
---|
| 294 | ! B.1. Get the ID of the variable |
---|
| 295 | ierr=NF_INQ_VARID(nid,name,varid) |
---|
| 296 | if (ierr.ne.NF_NOERR) then |
---|
| 297 | ! If we failed geting the variable's ID, we assume it is because |
---|
| 298 | ! the variable doesn't exist yet and must be created. |
---|
| 299 | ! Start by obtaining corresponding dimensions IDs |
---|
| 300 | ierr=NF_INQ_DIMID(nid,"time",id(1)) |
---|
| 301 | ! Tell the world about it |
---|
| 302 | write(*,*) "=====================" |
---|
| 303 | write(*,*) "writediagsoil: creating variable "//trim(name) |
---|
| 304 | call def_var(nid,name,title,units,1,id,varid,ierr) |
---|
| 305 | endif ! of if (ierr.ne.NF_NOERR) |
---|
| 306 | |
---|
| 307 | ! B.2. Prepare things to be able to write/append the variable |
---|
| 308 | corners(1)=ntime |
---|
| 309 | |
---|
| 310 | edges(1)=1 |
---|
| 311 | |
---|
| 312 | ! B.3. Write the data |
---|
| 313 | !#ifdef NC_DOUBLE |
---|
| 314 | ! ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0) |
---|
| 315 | !#else |
---|
| 316 | ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data0) |
---|
| 317 | !#endif |
---|
| 318 | if (ierr.ne.NF_NOERR) then |
---|
| 319 | write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//& |
---|
| 320 | " to file "//trim(filename)//" at time",date |
---|
| 321 | endif |
---|
| 322 | |
---|
| 323 | endif ! of if (dimpx.eq.3) elseif (dimpx.eq.2) ... |
---|
| 324 | endif ! of if (mod(zitau+1,isample).eq.0) |
---|
| 325 | |
---|
| 326 | ! 4. Close the NetCDF file |
---|
[1130] | 327 | if (is_master) then |
---|
| 328 | ierr=NF_CLOSE(nid) |
---|
| 329 | endif |
---|
[38] | 330 | |
---|
| 331 | end subroutine writediagsoil |
---|