source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 737

Last change on this file since 737 was 712, checked in by tnavarro, 13 years ago

bug if icelocationmode=1 + sanity check for ice negative values
File size: 18.2 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      IMPLICIT NONE
6c=======================================================================
7c
8c   creation des calottes pour l'etat initial
9c
10c=======================================================================
11c-----------------------------------------------------------------------
12c   Declarations:
13c   -------------
14#include "dimensions.h"
15#include "dimphys.h"
16#include "surfdat.h"
17#include "callkeys.h"
18#include "tracer.h"
19#include "comgeomfi.h"
20#include "comcstfi.h"
21
22#include "datafile.h"
23
24      INTEGER ngrid,ig,icap,iq,alternate
25      REAL  piceco2(ngrid),psolaralb(ngrid,2)
26      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
27      REAL icedryness ! ice dryness
28     
29      ! longwatercaptag is watercaptag. Trick for some compilers
30      LOGICAL, DIMENSION(100000) :: longwatercaptag
31
32      EXTERNAL ISMIN,ISMAX
33      INTEGER ISMIN,ISMAX
34     
35! There are 3 different modes for ice distribution:
36! icelocationmode = 1 ---> based on data from surface.nc
37! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
38! icelocationmode = 3 ---> based on logical relations for latitude and longitude
39! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
40      INTEGER,SAVE :: icelocationmode = 2
41       
42       
43      !in case icelocationmode == 1
44      INTEGER i,j
45      INTEGER     imd,jmd
46      PARAMETER   (imd=360,jmd=180)
47      REAL        zdata(imd,jmd)
48      REAL        zelat,zelon
49
50      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
51      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
52
53      REAL step,count,ratiolat
54
55      INTEGER   ierr,nid,nvarid
56     
57      REAL,SAVE :: min_icevalue = 500.
58      character(len=50) :: string = 'thermal'
59     
60      character (len=100) :: zedatafile
61c
62c=======================================================================
63
64c     water ice outliers
65c     ------------------------------------------
66
67      IF ((water) .and. (caps)) THEN
68     
69c Perennial H20 north cap defined by watercaptag=true (allows surface to be
70c hollowed by sublimation in vdifc).
71
72c We might not want albedodat to be modified because it is used to write
73c restart files. Instead, albedo is directly modified when needed (i.e.
74c if we have watercaptag and no co2 ice), below and in albedocaps.F90
75
76c       "Dryness coefficient" controlling the evaporation and
77c        sublimation from the ground water ice (close to 1)
78c        HERE, the goal is to correct for the fact
79c        that the simulated permanent water ice polar caps
80c        is larger than the actual cap and the atmospheric
81c        opacity not always realistic.
82
83         alternate = 0
84         
85         if (ngridmx .ne. 1) then
86           watercaptag(:) = .false.
87           longwatercaptag(:) = .false.
88         endif
89         
90         write(*,*) "Ice dryness ?"
91         icedryness=1. ! default value
92         call getin("icedryness",icedryness)
93         write(*,*) " icedryness = ",icedryness
94         dryness (:) = icedryness
95         
96       
97#ifdef MESOSCALE
98
99      do ig=1,ngridmx
100
101         !write(*,*) "all qsurf to zero. dirty."
102         do iq=1,nqmx
103         qsurf(ig,iq)=0.  !! on jette les inputs GCM
104                          !! on regle juste watercaptag
105                          !! il faudrait garder les inputs GCM
106                          !! si elles sont consequentes
107         enddo
108         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
109     .        ( albedodat(ig) .ge. 0.26   ) )  then
110                 write(*,*)"outlier ",ig
111                 watercaptag(ig)  = .true.
112                 dryness(ig)      = 1.
113                 albedodat(ig)    = albedo_h2o_ice  !! pour output
114         else
115                 watercaptag(ig)  = .false.
116                 dryness(ig)      = 1.
117         endif
118         
119      enddo
120#else
121
122
123
124      IF (ngridmx .eq. 1) THEN ! special case for 1d --> do nothing
125     
126         print*, 'ngridmx = 1, do no put ice caps in surfini.F'
127
128      ELSE IF (icelocationmode .eq. 1) THEN
129     
130         print*,'Surfini: ice caps defined from surface.nc'
131           
132! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
133! Typically, it is for thermal inertia above 500 tiu.
134! Two conditions are verified:
135! 1. GCM ice caps are defined such as area is conserved for a given latitude
136! (the approximation is that all points within the GCM latitude resolution have the same area).
137! 2. caps are placed to fill the GCM points with the most detected ice first.
138     
139
140           
141        zedatafile = trim(datafile)
142 
143       
144        ierr=nf90_open(trim(zedatafile)//'/surface.nc',
145     &  NF90_NOWRITE,nid)
146     
147      IF (ierr.NE.nf90_noerr) THEN
148       write(*,*)'Error : cannot open file surface.nc '
149       write(*,*)'(in phymars/surfini.F)'
150       write(*,*)'It should be in :',trim(zedatafile),'/'
151       write(*,*)'1) You can set this path in the callphys.def file:'
152       write(*,*)'   datadir=/path/to/the/datafiles'
153       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
154       write(*,*)'   can be obtained online on:'
155       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
156       CALL ABORT
157      ENDIF
158     
159     
160      ierr=nf90_inq_varid(nid, string, nvarid)
161      if (ierr.ne.nf90_noerr) then
162        write(*,*) 'surfini error, cannot find ',trim(string)
163        write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
164        write(*,*)trim(nf90_strerror(ierr))
165        stop
166      endif
167
168      ierr=nf90_get_var(nid, nvarid, zdata)
169
170      if (ierr.ne.nf90_noerr) then
171        write(*,*) 'surfini: error failed loading ',trim(string)
172        write(*,*)trim(nf90_strerror(ierr))
173        stop
174      endif
175 
176                     
177      ierr=nf90_close(nid)
178 
179
180      nb_ice(:,1) = 1 ! default: there is no ice
181      latice(:,1) = 1
182      lonice(:,1) = 1
183      nb_ice(:,2) = 0
184      latice(:,2) = 0
185      lonice(:,2) = 0
186      !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
187
188      ! loop over the GCM grid - except for poles (ig=1 and ngridmx)
189      do ig=2,ngridmx-1
190     
191        ! loop over the surface file grid     
192        do i=1,imd
193          do j=1,jmd
194            zelon = i - 180.
195            zelat = 90. - j
196            if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
197     &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
198              ! count all points in that GCM grid point
199              nb_ice(ig,1) = nb_ice(ig,1) + 1
200              if (zdata(i,j) > min_icevalue)
201                 ! count all detected points in that GCM grid point
202     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
203             endif
204          enddo
205        enddo 
206
207        ! projection of nb_ice on GCM lat and lon axes
208        latice(1+(ig-2)/iim,:) =
209     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
210        lonice(1+mod(ig-2,iim),:) =
211     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
212
213      enddo ! of do ig=2,ngridmx-1
214     
215
216     
217      ! special case for poles
218      nb_ice(1,2)   = 1  ! ice prescribed on north pole
219      latice(1,:)   = nb_ice(1,:)
220      lonice(1,:)   = nb_ice(1,:)
221      latice(jjm,:) = nb_ice(ngridmx,:)
222      lonice(iim,:) = nb_ice(ngridmx,:)
223     
224     
225!      print*, 'latice TOT', latice(:,1)
226!      print*, 'latice FOUND', latice(:,2)
227!      print*, 'lonice TOT', lonice(:,1)
228!      print*, 'lonice FOUND', lonice(:,2)
229     
230!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
231!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
232     
233!      print*,''
234!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
235!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
236     
237   
238      ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
239      do i=1,jjm/2
240      step  = 1. ! threshold to add ice cap
241      count = 0. ! number of ice GCM caps at this latitude
242      ! ratiolat is the ratio of area covered by ice within this GCM latitude range
243      ratiolat  = real(latice(i,2))/real(latice(i,1))
244      !print*,'i',i,(i-1)*iim+2,i*iim+1
245     
246        ! put ice caps while there is not enough ice,
247        ! as long as the threshold is above 20%
248        do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
249          count = 0.
250          ! loop over GCM longitudes
251          do j=1,iim
252            ! if the detected ice ratio in the GCM grid point
253            ! is more than 'step', then add ice
254            if (real(nb_ice((i-1)*iim+1+j,2))
255     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
256                  watercaptag((i-1)*iim+1+j) = .true.
257                  count = count + 1
258            endif
259          enddo ! of do j=1,iim
260          !print*, 'step',step,count,ratiolat*iim
261          step = step - 0.01
262        enddo ! of do while
263      !print*, 'step',step,count,ratiolat*iim
264
265      enddo ! of do i=1,jjm/2
266           
267
268      ELSE IF (icelocationmode .eq. 2) THEN
269     
270        print*,'Surfini: predefined ice caps'
271     
272        if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
273           
274          print*,'water ice caps distribution for 32x24 resolution'
275          longwatercaptag(1:9)    = .true. ! central cap - core
276          longwatercaptag(26:33)  = .true. ! central cap
277          longwatercaptag(1:33)  = .true. ! central cap
278          longwatercaptag(56)  = .true. ! central cap
279          longwatercaptag(58)  = .true. ! central cap
280          longwatercaptag(60)  = .true. ! central cap
281          longwatercaptag(62)  = .true. ! central cap
282          longwatercaptag(64)  = .true. ! central cap
283!---------------------   OUTLIERS  ----------------------------
284
285        else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
286
287          print*,'water ice caps distribution for 64x48 resolution'
288          longwatercaptag(1:65)   = .true. ! central cap - core
289          longwatercaptag(75:85)  = .true. ! central cap
290          longwatercaptag(93:114) = .true. ! central cap
291!---------------------   OUTLIERS  ----------------------------
292          if (.true.) then
293          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
294          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
295          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
296          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
297          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
298          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
299          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
300          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
301          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
302          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
303          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
304          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
305          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
306          longwatercaptag(194)    = .true. ! outlier, lat = 75
307          longwatercaptag(203)    = .true. ! outlier, lat = 75
308          longwatercaptag(207)    = .true. ! outlier, lat = 75
309          longwatercaptag(242)    = .true. ! outlier, lat = 75
310          longwatercaptag(244)    = .true. ! outlier, lat = 75
311          longwatercaptag(246)    = .true. ! outlier, lat = 75
312          longwatercaptag(248)    = .true. ! outlier, lat = 75
313          longwatercaptag(250)    = .true. ! outlier, lat = 75
314          longwatercaptag(252)    = .true. ! outlier, lat = 75
315          longwatercaptag(254)    = .true. ! outlier, lat = 75
316          longwatercaptag(256:257)= .true. ! outlier, lat = 75
317          endif
318!--------------------------------------------------------------       
319
320
321           
322        else if (ngridmx .ne. 1) then
323       
324      print*,'No predefined ice location for this resolution :',iim,jjm
325      print*,'Please change icelocationmode in surfini.F'
326      print*,'Or add some new definitions ...'
327      call abort
328         
329        endif
330
331        do ig=1,ngridmx
332          if (longwatercaptag(ig)) watercaptag(ig) = .true.
333        enddo
334
335
336      ELSE IF (icelocationmode .eq. 3) THEN
337     
338        print*,'Surfini: ice caps defined by lat and lon values'
339
340         do ig=1,ngridmx
341         
342c-------- Towards olympia planitia water caps -----------
343c--------------------------------------------------------
344
345       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
346     .           ( lati(ig)*180./pi .le. 80.  ) .and.
347     .           ( long(ig)*180./pi .ge. 110. ) .and.
348     .           ( long(ig)*180./pi .le. 181. ) )
349     .         .or.
350
351     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
352     .           ( lati(ig)*180./pi .le. 76.  ) .and.
353     .           ( long(ig)*180./pi .ge. 150. ) .and.
354     .           ( long(ig)*180./pi .le. 168. ) )
355     .         .or.
356     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
357     .           ( lati(ig)*180./pi .le. 80.  ) .and.
358     .           ( long(ig)*180./pi .ge. -150.) .and.
359     .           ( long(ig)*180./pi .le. -110.) ) )
360     .         then
361             
362               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
363              !    watercaptag(ig)=.true.
364                  alternate = 1
365               else
366                  alternate = 0
367               endif !end if alternate = 0
368               
369       endif
370
371c----------- Opposite olympia planitia water cap --------
372c--------------------------------------------------------
373
374        if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
375     .         ( lati(ig)*180./pi     .le.  84 ) )
376     .         .and.
377     .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
378     .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
379!     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
380!     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
381!     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
382!     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
383        !   watercaptag(ig)=.true.
384        endif
385
386
387c -------------------- Central cap ----------------------
388c--------------------------------------------------------
389
390      if (abs(lati(ig)*180./pi).gt.80)
391     .       watercaptag(ig)=.true.
392           
393c--------------------------------------------------------
394c--------------------------------------------------------
395      end do ! of (ngridmx)
396
397
398       ELSE
399     
400         print*, 'In surfini.F, icelocationmode is ', icelocationmode
401         print*, 'It should be 1, 2 or 3.'
402         call abort
403
404       ENDIF ! of if (icelocation)
405       
406       
407        ! print caps locations - useful for plots too
408        print*,'latitude | longitude | ig'
409        do ig=1,ngridmx
410          dryness (ig) = icedryness
411
412          if (watercaptag(ig)) then
413             print*,'ice water cap', lati(ig)*180./pi,
414     .              long(ig)*180./pi, ig
415          endif
416        enddo
417       
418#endif     
419
420       ENDIF ! (caps & water)
421       
422
423c ===============================================================
424c      INITIAL ALBEDO
425c ===============================================================
426
427         write(*,*)"surfini: water frost thickness",
428     s     frost_albedo_threshold
429         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
430         write(*,*)"surfini: water ice TI:", inert_h2o_ice
431
432c        To start with : Initial albedo = observed dataset
433c        -------------------------------------------------
434         DO ig=1,ngrid
435              psolaralb(ig,1)=albedodat(ig)
436              psolaralb(ig,2)=albedodat(ig)
437         END DO
438         PRINT*,'minimum albedo sans water caps',
439     s     albedodat(ISMIN(ngrid,albedodat,1))
440         PRINT*,'maximum albedo sans water caps',
441     s     albedodat(ISMAX(ngrid,albedodat,1))
442
443c        initial albedo if permanent H2O ice is present
444c        ------------------------------------------------
445         IF ((water) .and. (caps)) THEN
446           DO ig=1,ngrid
447            IF (watercaptag(ig)) THEN
448              psolaralb(ig,1) = albedo_h2o_ice
449              psolaralb(ig,2) = albedo_h2o_ice
450            ENDIF
451           END DO
452           PRINT*,'minimum albedo avec water caps',
453     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
454           PRINT*,'maximum albedo avec water caps',
455     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
456
457         ENDIF
458
459c      changing initial albedo if CO2 ice is present
460c      -------------------------------------------
461
462       DO ig=1,ngrid
463         IF (piceco2(ig) .GT. 0.) THEN
464             IF(ig.GT.ngrid/2+1) THEN
465                icap=2
466             ELSE
467                icap=1
468             ENDIF
469             psolaralb(ig,1) = albedice(icap)
470             psolaralb(ig,2) = albedice(icap)   
471         END IF
472       END DO
473
474c      changing initial albedo if water ice frost is present
475c      -------------------------------------------
476       IF (water) THEN
477          do iq=1,nqmx
478c          if there is frost and surface albedo is set to albedo_h2o_ice
479           if(noms(iq).eq."h2o_ice") then
480             do ig=1,ngrid
481             
482              if ((watercaptag(ig).eqv..false.)
483     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
484              print*, ''
485              print*, '!!! PROBLEM in SURFINI !!!!'
486              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
487     & WATERCAPTAG IS FALSE'
488              print*, ''
489              print*, 'ig,qsurf,threshold' ,
490     &         ig, qsurf(ig,iq), -frost_albedo_threshold
491              print*, ''
492              print*, '1) Check h2o_ice in startfi and ice
493     & distribution in surfini'
494              print*, '2) Use ini_h2osurf option in newstart'
495              print*, ''
496              CALL ABORT
497            endif
498             
499              if ((piceco2(ig) .eq. 0.).and.
500     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
501                     psolaralb(ig,1) = albedo_h2o_ice
502                     psolaralb(ig,2) = albedo_h2o_ice
503c                     PRINT*,'surfini.F frost',
504c     &                  lati(ig)*180./pi, long(ig)*180./pi
505               endif
506              enddo
507           endif
508          end do
509          PRINT*,'minimum albedo avec givre et co2',
510     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
511          PRINT*,'maximum albedo avec givre et co2',
512     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
513       END IF
514         
515      RETURN
516      END
Note: See TracBrowser for help on using the repository browser.