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

Last change on this file since 1087 was 1087, checked in by tnavarro, 11 years ago

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